qualia_core_db/solvers/statistics/
information.rs

1//! Information theory — Shannon entropy, KL divergence, cross-entropy and mutual
2//! information over discrete distributions / samples (all in **bits**, `log₂`).
3//!
4//! Mission note: mutual information is a principled, assumption-free relevance
5//! signal — `I(X;Y)` measures how much knowing `X` reduces uncertainty about `Y`
6//! with no linearity assumption — which is exactly what the 10D→5D NQuin relevance
7//! router needs to choose its projection.
8
9/// Shannon entropy `H(p) = −Σ pᵢ·log₂ pᵢ` (bits) of a probability vector. Zero
10/// probabilities contribute 0. `None` if empty or the masses don't form a positive
11/// distribution.
12pub fn entropy(p: &[f64]) -> Option<f64> {
13    if p.is_empty() {
14        return None;
15    }
16    let total: f64 = p.iter().sum();
17    if !(total > 0.0) {
18        return None;
19    }
20    let mut h = 0.0;
21    for &pi in p {
22        let q = pi / total; // tolerate unnormalized input
23        if q > 0.0 {
24            h -= q * q.log2();
25        }
26    }
27    Some(h)
28}
29
30/// Entropy from integer counts (normalized internally).
31pub fn entropy_from_counts(counts: &[usize]) -> Option<f64> {
32    if counts.is_empty() {
33        return None;
34    }
35    let p: Vec<f64> = counts.iter().map(|&c| c as f64).collect();
36    entropy(&p)
37}
38
39/// Kullback–Leibler divergence `D(p‖q) = Σ pᵢ·log₂(pᵢ/qᵢ)` (bits). Both inputs are
40/// normalized internally. `None` on a length mismatch, empty input, or if `qᵢ = 0`
41/// where `pᵢ > 0` (the divergence is then infinite — refuse rather than fabricate).
42pub fn kl_divergence(p: &[f64], q: &[f64]) -> Option<f64> {
43    if p.is_empty() || p.len() != q.len() {
44        return None;
45    }
46    let (sp, sq): (f64, f64) = (p.iter().sum(), q.iter().sum());
47    if !(sp > 0.0) || !(sq > 0.0) {
48        return None;
49    }
50    let mut d = 0.0;
51    for (&pi, &qi) in p.iter().zip(q) {
52        let pn = pi / sp;
53        let qn = qi / sq;
54        if pn > 0.0 {
55            if qn <= 0.0 {
56                return None; // support of p not covered by q
57            }
58            d += pn * (pn / qn).log2();
59        }
60    }
61    Some(d)
62}
63
64/// Cross-entropy `H(p, q) = −Σ pᵢ·log₂ qᵢ` (bits). `None` like [`kl_divergence`].
65pub fn cross_entropy(p: &[f64], q: &[f64]) -> Option<f64> {
66    Some(entropy(p)? + kl_divergence(p, q)?)
67}
68
69/// Mutual information `I(X;Y) = H(X) + H(Y) − H(X,Y)` (bits), estimated from paired
70/// discrete samples (small non-negative integer labels). `None` on length mismatch
71/// or empty input. `I ≥ 0`, and `I = 0` iff `X ⟂ Y` in the sample.
72pub fn mutual_information_discrete(x: &[usize], y: &[usize]) -> Option<f64> {
73    let n = x.len();
74    if n == 0 || n != y.len() {
75        return None;
76    }
77    let nx = x.iter().max().copied().unwrap_or(0) + 1;
78    let ny = y.iter().max().copied().unwrap_or(0) + 1;
79    let mut joint = vec![0.0f64; nx * ny];
80    let mut px = vec![0.0f64; nx];
81    let mut py = vec![0.0f64; ny];
82    for (&xi, &yi) in x.iter().zip(y) {
83        joint[xi * ny + yi] += 1.0;
84        px[xi] += 1.0;
85        py[yi] += 1.0;
86    }
87    let nf = n as f64;
88    let mut mi = 0.0;
89    for xi in 0..nx {
90        for yi in 0..ny {
91            let pxy = joint[xi * ny + yi] / nf;
92            if pxy > 0.0 {
93                let pxi = px[xi] / nf;
94                let pyi = py[yi] / nf;
95                mi += pxy * (pxy / (pxi * pyi)).log2();
96            }
97        }
98    }
99    Some(mi.max(0.0))
100}
101
102#[cfg(test)]
103mod tests {
104    use super::*;
105    const EPS: f64 = 1e-9;
106
107    #[test]
108    fn entropy_known_values() {
109        // Fair coin → 1 bit; fair 4-way → 2 bits; certain outcome → 0.
110        assert!((entropy(&[0.5, 0.5]).unwrap() - 1.0).abs() < EPS);
111        assert!((entropy(&[0.25; 4]).unwrap() - 2.0).abs() < EPS);
112        assert!(entropy(&[1.0, 0.0, 0.0]).unwrap().abs() < EPS);
113        // Unnormalized counts work too.
114        assert!((entropy_from_counts(&[1, 1, 1, 1]).unwrap() - 2.0).abs() < EPS);
115    }
116
117    #[test]
118    fn kl_is_zero_for_equal_and_positive_otherwise() {
119        assert!(kl_divergence(&[0.5, 0.5], &[0.5, 0.5]).unwrap().abs() < EPS);
120        assert!(kl_divergence(&[0.9, 0.1], &[0.5, 0.5]).unwrap() > 0.0);
121        // Infinite divergence (q has no support where p does) → refuse.
122        assert!(kl_divergence(&[0.5, 0.5], &[1.0, 0.0]).is_none());
123    }
124
125    #[test]
126    fn mutual_information_detects_dependence() {
127        // y = x → I(X;Y) = H(X) = 1 bit for a balanced binary x.
128        let x = [0usize, 0, 1, 1, 0, 1, 0, 1];
129        let y = x;
130        assert!((mutual_information_discrete(&x, &y).unwrap() - 1.0).abs() < EPS);
131        // Independent x,y → MI ≈ 0.
132        let xi = [0usize, 0, 1, 1, 0, 0, 1, 1];
133        let yi = [0usize, 1, 0, 1, 0, 1, 0, 1];
134        assert!(mutual_information_discrete(&xi, &yi).unwrap() < 1e-9);
135    }
136
137    #[test]
138    fn cross_entropy_decomposes() {
139        // H(p,q) = H(p) + D(p||q).
140        let p = [0.7, 0.3];
141        let q = [0.5, 0.5];
142        let ce = cross_entropy(&p, &q).unwrap();
143        let h = entropy(&p).unwrap();
144        let d = kl_divergence(&p, &q).unwrap();
145        assert!((ce - (h + d)).abs() < EPS);
146    }
147
148    #[test]
149    fn guards() {
150        assert_eq!(entropy(&[]), None);
151        assert_eq!(kl_divergence(&[0.5, 0.5], &[1.0]), None);
152        assert_eq!(mutual_information_discrete(&[0, 1], &[0]), None);
153    }
154}
qualia_core_db/solvers/statistics/information.rs

qualia_core_db/solvers/statistics/
information.rs
