sc_neurocore_engine/analysis/
patterns.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 directionality, ordering, and higher-order patterns
8
9/// Spike directionality (Kreuz et al. 2015).
10/// Returns asymmetry in [-1, 1]. Positive: A leads B.
11use rayon::prelude::*;
12
13pub fn spike_directionality(times_a: &[f64], times_b: &[f64], t_start: f64, t_end: f64) -> f64 {
14    let mut ta: Vec<f64> = times_a
15        .iter()
16        .copied()
17        .filter(|&t| t >= t_start && t <= t_end)
18        .collect();
19    let tb: Vec<f64> = times_b
20        .iter()
21        .copied()
22        .filter(|&t| t >= t_start && t <= t_end)
23        .collect();
24    ta.sort_by(|a, b| a.partial_cmp(b).unwrap());
25
26    if ta.is_empty() || tb.is_empty() {
27        return 0.0;
28    }
29
30    let mut lead_a = 0_usize;
31    let mut lead_b = 0_usize;
32
33    // Ensure tb is sorted for binary search
34    let mut tb_sorted = tb.clone();
35    tb_sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
36
37    for &t in &ta {
38        let idx = tb_sorted.partition_point(|&x| x < t);
39
40        let nearest_after = if idx < tb_sorted.len() {
41            Some(tb_sorted[idx] - t)
42        } else {
43            None
44        };
45
46        let nearest_before = if idx > 0 {
47            Some(t - tb_sorted[idx - 1])
48        } else {
49            None
50        };
51
52        if let (Some(nb), Some(na)) = (nearest_before, nearest_after) {
53            if nb < na {
54                lead_b += 1;
55            } else {
56                lead_a += 1;
57            }
58        }
59    }
60
61    let total = lead_a + lead_b;
62    if total == 0 {
63        return 0.0;
64    }
65    (lead_a as f64 - lead_b as f64) / total as f64
66}
67
68/// Spike train order matrix (Kreuz et al. 2017).
69/// Returns n×n matrix (flat, row-major) of pairwise directionality.
70pub fn spike_train_order(times_list: &[&[f64]], t_start: f64, t_end: f64) -> Vec<f64> {
71    let n = times_list.len();
72    let mut mat = vec![0.0_f64; n * n];
73    mat.par_chunks_exact_mut(n)
74        .enumerate()
75        .for_each(|(i, row)| {
76            for j in 0..n {
77                if i == j {
78                    continue;
79                }
80                // Note: full matrix computation simplifies parallel dispatch
81                // even though directionality is antisymmetric.
82                if j > i {
83                    let d = spike_directionality(times_list[i], times_list[j], t_start, t_end);
84                    row[j] = d;
85                } else {
86                    let d = spike_directionality(times_list[j], times_list[i], t_start, t_end);
87                    row[j] = -d;
88                }
89            }
90        });
91    mat
92}
93
94/// Third-order cumulant C3(tau1, tau2) (Nikias & Petropulu 1993).
95/// Returns max_lag × max_lag matrix (flat, row-major).
96pub fn cubic_higher_order(binary_train: &[i32], dt: f64, max_lag: usize) -> Vec<f64> {
97    let _ = dt;
98    let n = binary_train.len();
99    let mean: f64 = binary_train.iter().map(|&v| v as f64).sum::<f64>() / n as f64;
100    let x: Vec<f64> = binary_train.iter().map(|&v| v as f64 - mean).collect();
101
102    let mut c3 = vec![0.0_f64; max_lag * max_lag];
103    c3.par_chunks_exact_mut(max_lag)
104        .enumerate()
105        .for_each(|(t1, row)| {
106            for t2 in 0..max_lag {
107                let valid_n = n.saturating_sub(t1.max(t2));
108                if valid_n == 0 {
109                    continue;
110                }
111                let mut sum = 0.0_f64;
112                let mut k = 0;
113                while k + 3 < valid_n {
114                    sum += x[k] * x[k + t1] * x[k + t2];
115                    sum += x[k + 1] * x[k + 1 + t1] * x[k + 1 + t2];
116                    sum += x[k + 2] * x[k + 2 + t1] * x[k + 2 + t2];
117                    sum += x[k + 3] * x[k + 3 + t1] * x[k + 3 + t2];
118                    k += 4;
119                }
120                while k < valid_n {
121                    sum += x[k] * x[k + t1] * x[k + t2];
122                    k += 1;
123                }
124                row[t2] = sum / valid_n as f64;
125            }
126        });
127    c3
128}
129
130#[cfg(test)]
131mod tests {
132    use super::*;
133
134    // ── spike_directionality ────────────────────────────────────────
135
136    #[test]
137    fn test_directionality_a_leads() {
138        // A fires before B consistently
139        let ta = vec![0.1, 0.3, 0.5, 0.7];
140        let tb = vec![0.15, 0.35, 0.55, 0.75];
141        let d = spike_directionality(&ta, &tb, 0.0, 1.0);
142        assert!(d > 0.0, "A leads B → positive, got {d}");
143    }
144
145    #[test]
146    fn test_directionality_b_leads() {
147        let ta = vec![0.15, 0.35, 0.55, 0.75];
148        let tb = vec![0.1, 0.3, 0.5, 0.7];
149        let d = spike_directionality(&ta, &tb, 0.0, 1.0);
150        assert!(d < 0.0, "B leads A → negative, got {d}");
151    }
152
153    #[test]
154    fn test_directionality_empty() {
155        assert_eq!(spike_directionality(&[], &[0.5], 0.0, 1.0), 0.0);
156    }
157
158    #[test]
159    fn test_directionality_antisymmetric() {
160        let ta = vec![0.1, 0.3, 0.5];
161        let tb = vec![0.15, 0.35, 0.55];
162        let d_ab = spike_directionality(&ta, &tb, 0.0, 1.0);
163        let d_ba = spike_directionality(&tb, &ta, 0.0, 1.0);
164        // Not exactly antisymmetric due to algorithm, but should have opposite signs
165        assert!(d_ab > 0.0 && d_ba < 0.0, "should have opposite signs");
166    }
167
168    // ── spike_train_order ───────────────────────────────────────────
169
170    #[test]
171    fn test_order_antisymmetric() {
172        let ta = vec![0.1, 0.3, 0.5];
173        let tb = vec![0.15, 0.35, 0.55];
174        let trains: Vec<&[f64]> = vec![&ta, &tb];
175        let mat = spike_train_order(&trains, 0.0, 1.0);
176        assert!((mat[1] + mat[2]).abs() < 1e-10, "antisymmetric");
177        assert_eq!(mat[0], 0.0, "diagonal = 0");
178    }
179
180    #[test]
181    fn test_order_shape() {
182        let ta = vec![0.1];
183        let tb = vec![0.2];
184        let tc = vec![0.3];
185        let trains: Vec<&[f64]> = vec![&ta, &tb, &tc];
186        let mat = spike_train_order(&trains, 0.0, 1.0);
187        assert_eq!(mat.len(), 9);
188    }
189
190    // ── cubic_higher_order ──────────────────────────────────────────
191
192    #[test]
193    fn test_c3_shape() {
194        let train = vec![0i32, 1, 0, 1, 0, 0, 1, 0, 1, 0];
195        let c3 = cubic_higher_order(&train, 0.001, 5);
196        assert_eq!(c3.len(), 25);
197    }
198
199    #[test]
200    fn test_c3_zero_lag_positive() {
201        // C3(0,0) = E[x^3] which is positive for sparse spikes (skewed)
202        let train = vec![0i32, 0, 0, 0, 0, 0, 0, 0, 0, 1];
203        let c3 = cubic_higher_order(&train, 0.001, 3);
204        assert!(
205            c3[0] > 0.0,
206            "C3(0,0) should be positive for skewed data, got {}",
207            c3[0]
208        );
209    }
210
211    #[test]
212    fn test_c3_constant_zero() {
213        let train = vec![0i32; 100];
214        let c3 = cubic_higher_order(&train, 0.001, 5);
215        assert!(c3.iter().all(|&v| v.abs() < 1e-10), "constant → C3 = 0");
216    }
217}
sc_neurocore_engine/analysis/patterns.rs

sc_neurocore_engine/analysis/
patterns.rs
