sc_neurocore_engine/
phi.rs

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

sc_neurocore_engine/
phi.rs
