sc_neurocore_engine/
phi.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 — Phi* integrated information estimation in Rust
8
9//! Barrett & Seth 2011 geometric Phi under Gaussian assumption.
10//! Phi* = MI(past; future) - max_partition sum MI(past_k; future_k)
11
12/// Compute Phi* for a multi-channel time series.
13/// data: [n_channels][n_timesteps], row-major.
14pub fn phi_star(data: &[Vec<f64>], tau: usize) -> f64 {
15    let n = data.len();
16    if n < 2 || data[0].len() <= 2 * tau {
17        return 0.0;
18    }
19    let t = data[0].len();
20
21    // Split into past and future
22    let past: Vec<Vec<f64>> = data.iter().map(|ch| ch[..t - tau].to_vec()).collect();
23    let future: Vec<Vec<f64>> = data.iter().map(|ch| ch[tau..].to_vec()).collect();
24
25    let mi_whole = gaussian_mi(&past, &future);
26
27    // Minimum information partition: contiguous splits only
28    let mut mi_parts_min = f64::INFINITY;
29    for k in 1..n {
30        let past_a: Vec<Vec<f64>> = past[..k].to_vec();
31        let fut_a: Vec<Vec<f64>> = future[..k].to_vec();
32        let past_b: Vec<Vec<f64>> = past[k..].to_vec();
33        let fut_b: Vec<Vec<f64>> = future[k..].to_vec();
34        let mi_parts = gaussian_mi(&past_a, &fut_a) + gaussian_mi(&past_b, &fut_b);
35        if mi_parts < mi_parts_min {
36            mi_parts_min = mi_parts;
37        }
38    }
39
40    (mi_whole - mi_parts_min).max(0.0)
41}
42
43/// Gaussian mutual information: MI(X;Y) = 0.5 * ln(det(Cov_X) * det(Cov_Y) / det(Cov_XY))
44fn gaussian_mi(x: &[Vec<f64>], y: &[Vec<f64>]) -> f64 {
45    let eps = 1e-10;
46    let cov_x = covariance_matrix(x, eps);
47    let cov_y = covariance_matrix(y, eps);
48
49    let mut xy: Vec<Vec<f64>> = x.to_vec();
50    xy.extend_from_slice(y);
51    let cov_xy = covariance_matrix(&xy, eps);
52
53    let det_x = determinant(&cov_x).max(1e-300);
54    let det_y = determinant(&cov_y).max(1e-300);
55    let det_xy = determinant(&cov_xy).max(1e-300);
56
57    let mi = 0.5 * (det_x * det_y / det_xy).ln();
58    mi.max(0.0)
59}
60
61/// Compute covariance matrix with regularization.
62fn covariance_matrix(data: &[Vec<f64>], eps: f64) -> Vec<Vec<f64>> {
63    let n = data.len();
64    if n == 0 {
65        return vec![];
66    }
67    let t = data[0].len();
68    if t == 0 {
69        return vec![vec![eps; n]; n];
70    }
71
72    // Compute means
73    let means: Vec<f64> = data
74        .iter()
75        .map(|ch| ch.iter().sum::<f64>() / t as f64)
76        .collect();
77
78    // Compute covariance
79    let mut cov = vec![vec![0.0; n]; n];
80    let t_f = t as f64;
81    let denom = (t_f - 1.0).max(1.0);
82
83    for i in 0..n {
84        for j in 0..=i {
85            let dot = crate::simd::dot_f64_dispatch(&data[i], &data[j]);
86            let val = (dot - t_f * means[i] * means[j]) / denom;
87            cov[i][j] = val;
88            cov[j][i] = val;
89        }
90        cov[i][i] += eps; // regularize diagonal
91    }
92    cov
93}
94
95/// Determinant via LU decomposition (in-place, partial pivoting).
96fn determinant(matrix: &[Vec<f64>]) -> f64 {
97    let n = matrix.len();
98    if n == 0 {
99        return 1.0;
100    }
101    if n == 1 {
102        return matrix[0][0];
103    }
104    if n == 2 {
105        return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
106    }
107
108    let mut a: Vec<Vec<f64>> = matrix.to_vec();
109    let mut det = 1.0f64;
110
111    for col in 0..n {
112        // Partial pivot
113        let mut max_row = col;
114        let mut max_val = a[col][col].abs();
115        for row in (col + 1)..n {
116            if a[row][col].abs() > max_val {
117                max_val = a[row][col].abs();
118                max_row = row;
119            }
120        }
121        if max_val < 1e-300 {
122            return 0.0;
123        }
124        if max_row != col {
125            a.swap(col, max_row);
126            det = -det;
127        }
128
129        det *= a[col][col];
130        let pivot = a[col][col];
131        for row in (col + 1)..n {
132            let factor = a[row][col] / pivot;
133            for j in (col + 1)..n {
134                a[row][j] -= factor * a[col][j];
135            }
136        }
137    }
138    det
139}
140
141#[cfg(test)]
142mod tests {
143    use super::*;
144
145    #[test]
146    fn test_independent_channels_low_phi() {
147        // Independent random channels: Phi should be near zero
148        let data = vec![
149            (0..200).map(|i| (i as f64 * 0.1).sin()).collect(),
150            (0..200).map(|i| (i as f64 * 0.3).cos()).collect(),
151            (0..200).map(|i| (i as f64 * 0.7).sin()).collect(),
152        ];
153        let phi = phi_star(&data, 1);
154        assert!(phi < 0.5, "Expected low Phi for independent, got {phi}");
155    }
156
157    #[test]
158    fn test_correlated_channels_positive_phi() {
159        // Cross-correlated but distinct: Phi should be positive
160        let data = vec![
161            (0..200)
162                .map(|i| (i as f64 * 0.1).sin() + (i as f64 * 0.05).cos() * 0.3)
163                .collect(),
164            (0..200)
165                .map(|i| (i as f64 * 0.1).sin() * 0.8 + (i as f64 * 0.2).sin() * 0.5)
166                .collect(),
167            (0..200)
168                .map(|i| (i as f64 * 0.1).sin() * 0.6 + (i as f64 * 0.15).cos() * 0.7)
169                .collect(),
170        ];
171        let phi = phi_star(&data, 1);
172        assert!(phi >= 0.0, "Expected non-negative Phi, got {phi}");
173    }
174
175    #[test]
176    fn test_single_channel_zero() {
177        let data = vec![vec![1.0, 2.0, 3.0, 4.0, 5.0]];
178        assert_eq!(phi_star(&data, 1), 0.0);
179    }
180
181    #[test]
182    fn test_short_data_zero() {
183        let data = vec![vec![1.0, 2.0], vec![3.0, 4.0]];
184        assert_eq!(phi_star(&data, 1), 0.0);
185    }
186
187    #[test]
188    fn test_determinant_2x2() {
189        let m = vec![vec![1.0, 2.0], vec![3.0, 4.0]];
190        let det = determinant(&m);
191        assert!((det - (-2.0)).abs() < 1e-10);
192    }
193
194    #[test]
195    fn test_determinant_3x3() {
196        let m = vec![
197            vec![1.0, 2.0, 3.0],
198            vec![0.0, 1.0, 4.0],
199            vec![5.0, 6.0, 0.0],
200        ];
201        let det = determinant(&m);
202        assert!((det - 1.0).abs() < 1e-10);
203    }
204}
sc_neurocore_engine/phi.rs

sc_neurocore_engine/
phi.rs
