sc_neurocore_engine/analysis/
waveform.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 — Spike waveform shape analysis
8
9/// Trough-to-peak width in seconds. Bartho et al. 2004.
10pub fn waveform_width(waveform: &[f64], dt: f64) -> f64 {
11    if waveform.is_empty() {
12        return f64::NAN;
13    }
14    let trough = waveform
15        .iter()
16        .enumerate()
17        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap())
18        .unwrap()
19        .0;
20    if trough >= waveform.len() - 1 {
21        return f64::NAN;
22    }
23    let peak = trough
24        + waveform[trough..]
25            .iter()
26            .enumerate()
27            .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
28            .unwrap()
29            .0;
30    (peak - trough) as f64 * dt
31}
32
33/// Peak-to-trough amplitude. Bartho et al. 2004.
34pub fn waveform_amplitude(waveform: &[f64]) -> f64 {
35    if waveform.is_empty() {
36        return 0.0;
37    }
38    let max = waveform.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
39    let min = waveform.iter().cloned().fold(f64::INFINITY, f64::min);
40    max - min
41}
42
43/// Repolarisation slope: max dV/dt after trough. Bean 2007.
44pub fn waveform_repolarization_slope(waveform: &[f64], dt: f64) -> f64 {
45    if waveform.is_empty() {
46        return f64::NAN;
47    }
48    let trough = waveform
49        .iter()
50        .enumerate()
51        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap())
52        .unwrap()
53        .0;
54    if trough >= waveform.len() - 2 {
55        return f64::NAN;
56    }
57    let post = &waveform[trough..];
58    let mut max_dv = f64::NEG_INFINITY;
59    for i in 0..post.len() - 1 {
60        let dv = (post[i + 1] - post[i]) / dt;
61        if dv > max_dv {
62            max_dv = dv;
63        }
64    }
65    max_dv
66}
67
68/// Recovery slope: dV/dt during return to baseline after peak. Bean 2007.
69pub fn waveform_recovery_slope(waveform: &[f64], dt: f64) -> f64 {
70    if waveform.is_empty() {
71        return f64::NAN;
72    }
73    let trough = waveform
74        .iter()
75        .enumerate()
76        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap())
77        .unwrap()
78        .0;
79    if trough >= waveform.len() - 1 {
80        return f64::NAN;
81    }
82    let peak = trough
83        + waveform[trough..]
84            .iter()
85            .enumerate()
86            .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
87            .unwrap()
88            .0;
89    if peak >= waveform.len() - 2 {
90        return f64::NAN;
91    }
92    let post = &waveform[peak..];
93    let mut min_dv = f64::INFINITY;
94    for i in 0..post.len() - 1 {
95        let dv = (post[i + 1] - post[i]) / dt;
96        if dv < min_dv {
97            min_dv = dv;
98        }
99    }
100    min_dv
101}
102
103/// Half-width: duration at half-minimum amplitude. Bartho et al. 2004.
104pub fn waveform_halfwidth(waveform: &[f64], dt: f64) -> f64 {
105    if waveform.is_empty() {
106        return f64::NAN;
107    }
108    let trough_val = waveform.iter().cloned().fold(f64::INFINITY, f64::min);
109    let half_val = trough_val / 2.0;
110    let below: Vec<usize> = waveform
111        .iter()
112        .enumerate()
113        .filter(|(_, &v)| v < half_val)
114        .map(|(i, _)| i)
115        .collect();
116    if below.len() < 2 {
117        return f64::NAN;
118    }
119    (below[below.len() - 1] - below[0]) as f64 * dt
120}
121
122/// Peak-to-trough ratio. Bartho et al. 2004.
123pub fn waveform_pt_ratio(waveform: &[f64]) -> f64 {
124    if waveform.is_empty() {
125        return f64::NAN;
126    }
127    let trough = waveform
128        .iter()
129        .enumerate()
130        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap())
131        .unwrap()
132        .0;
133    let trough_val = waveform[trough].abs();
134    if trough >= waveform.len() - 1 || trough_val < 1e-30 {
135        return f64::NAN;
136    }
137    let peak_val = waveform[trough..]
138        .iter()
139        .cloned()
140        .fold(f64::NEG_INFINITY, f64::max);
141    peak_val.abs() / trough_val
142}
143
144#[cfg(test)]
145mod tests {
146    use super::*;
147
148    fn typical_waveform() -> Vec<f64> {
149        // Simulated extracellular action potential shape
150        vec![
151            0.0, 0.1, 0.0, -0.3, -0.8, -1.0, -0.7, -0.2, 0.3, 0.5, 0.4, 0.2, 0.1, 0.0,
152        ]
153    }
154
155    #[test]
156    fn test_waveform_width() {
157        let wf = typical_waveform();
158        let dt = 1.0 / 30000.0;
159        let w = waveform_width(&wf, dt);
160        // Trough at index 5, peak after at index 9 -> 4 samples
161        assert!((w - 4.0 * dt).abs() < 1e-12);
162    }
163
164    #[test]
165    fn test_waveform_width_empty() {
166        assert!(waveform_width(&[], 1.0 / 30000.0).is_nan());
167    }
168
169    #[test]
170    fn test_waveform_amplitude() {
171        let wf = typical_waveform();
172        let amp = waveform_amplitude(&wf);
173        // max=0.5, min=-1.0 -> 1.5
174        assert!((amp - 1.5).abs() < 1e-12);
175    }
176
177    #[test]
178    fn test_waveform_amplitude_empty() {
179        assert_eq!(waveform_amplitude(&[]), 0.0);
180    }
181
182    #[test]
183    fn test_repolarization_slope() {
184        let wf = typical_waveform();
185        let dt = 1.0 / 30000.0;
186        let slope = waveform_repolarization_slope(&wf, dt);
187        assert!(slope > 0.0);
188        // Max positive derivative after trough at index 5
189        // Steepest rise: index 7->8: 0.3-(-0.2)=0.5
190        let expected = 0.5 / dt;
191        assert!((slope - expected).abs() < 1e-6);
192    }
193
194    #[test]
195    fn test_recovery_slope() {
196        let wf = typical_waveform();
197        let dt = 1.0 / 30000.0;
198        let slope = waveform_recovery_slope(&wf, dt);
199        // Should be negative (descending from peak)
200        assert!(slope < 0.0);
201    }
202
203    #[test]
204    fn test_halfwidth() {
205        let wf = typical_waveform();
206        let dt = 1.0 / 30000.0;
207        let hw = waveform_halfwidth(&wf, dt);
208        // trough_val = -1.0, half_val = -0.5
209        // Below -0.5: indices 4 (-0.8), 5 (-1.0), 6 (-0.7) -> span = (6-4)*dt = 2*dt
210        assert!((hw - 2.0 * dt).abs() < 1e-12);
211    }
212
213    #[test]
214    fn test_pt_ratio() {
215        let wf = typical_waveform();
216        let ratio = waveform_pt_ratio(&wf);
217        // trough_val = 1.0, peak after trough = 0.5
218        assert!((ratio - 0.5).abs() < 1e-12);
219    }
220
221    #[test]
222    fn test_pt_ratio_no_peak() {
223        let wf = vec![0.0, -1.0];
224        assert!(waveform_pt_ratio(&wf).is_nan());
225    }
226
227    #[test]
228    fn test_flat_waveform() {
229        let wf = vec![0.5; 10];
230        // flat -> min_by picks first index (0), max after it spans to end
231        // trough at 0, peak at 0 (all equal) -> width = 0
232        // Actually argmax of [0.5, 0.5, ...] picks idx 0 relative -> peak - trough = 0
233        let w = waveform_width(&wf, 1.0 / 30000.0);
234        // For flat waveform, trough=0, max from [0..] is also 0, so width = 0
235        assert!(w.abs() < 1e-6 || w >= 0.0, "Width for flat waveform: {w}");
236    }
237
238    #[test]
239    fn test_single_sample() {
240        assert!(waveform_width(&[1.0], 0.001).is_nan());
241        assert!((waveform_amplitude(&[5.0])).abs() < 1e-12);
242    }
243
244    #[test]
245    fn test_negative_only_waveform() {
246        let wf = vec![-0.5, -1.0, -0.3];
247        let dt = 1.0 / 30000.0;
248        let w = waveform_width(&wf, dt);
249        // trough at 1, peak after at 2 -> 1 sample
250        assert!((w - dt).abs() < 1e-12);
251    }
252}
sc_neurocore_engine/analysis/waveform.rs

sc_neurocore_engine/analysis/
waveform.rs
