flow_linalg/
hotspot.rs

1//! Mage Hotspot Matrix: \(H = \sqrt{|S^{-1}|}\) where \(S\) is cosine similarity.
2//!
3//! Diagonal entries are Spreading Inflation Factors (SIFs). Off-diagonal entries
4//! indicate fluorochrome combinations that drive or suffer unmixing-dependent spread.
5
6use faer::linalg::solvers::{DenseSolveCore, PartialPivLu};
7use faer::{Mat, MatRef};
8
9/// Hotspot matrix \(H_{ij} = \sqrt{|(S^{-1})_{ij}|}\).
10#[derive(Debug, Clone)]
11pub struct HotspotMatrix {
12    pub matrix: Mat<f64>,
13}
14
15impl HotspotMatrix {
16    /// Diagonal SIFs (one per endmember).
17    pub fn sifs(&self) -> Vec<f64> {
18        let n = self.matrix.nrows();
19        (0..n).map(|i| self.matrix[(i, i)]).collect()
20    }
21
22    /// Row-major \(n \times n\) flat buffer.
23    pub fn flat_row_major(&self) -> Vec<f64> {
24        let n = self.matrix.nrows();
25        let mut out = Vec::with_capacity(n * n);
26        for i in 0..n {
27            for j in 0..n {
28                out.push(self.matrix[(i, j)]);
29            }
30        }
31        out
32    }
33}
34
35/// Compute hotspot from a square cosine-similarity (Gram) matrix \(S\).
36pub fn hotspot_from_similarity(similarity: MatRef<'_, f64>) -> Result<HotspotMatrix, String> {
37    let n = similarity.nrows();
38    if n == 0 || similarity.ncols() != n {
39        return Err("hotspot requires a non-empty square similarity matrix".into());
40    }
41    let lu = PartialPivLu::new(similarity);
42    let u = lu.U();
43    for i in 0..n {
44        if !u[(i, i)].is_finite() || u[(i, i)].abs() < 1e-12 {
45            return Err(format!(
46                "similarity matrix is singular or ill-conditioned at diagonal index {i}"
47            ));
48        }
49    }
50    let inv = lu.inverse();
51    let matrix = Mat::<f64>::from_fn(n, n, |i, j| inv[(i, j)].abs().sqrt());
52    Ok(HotspotMatrix { matrix })
53}
54
55/// Unit-normalize mixing-matrix columns, form \(S = A_u^\top A_u\), then hotspot.
56pub fn hotspot_from_mixing_matrix(mixing: MatRef<'_, f64>) -> Result<HotspotMatrix, String> {
57    let m = mixing.nrows();
58    let n = mixing.ncols();
59    if m == 0 || n == 0 {
60        return Err("hotspot requires a non-empty mixing matrix".into());
61    }
62    let mut au = Mat::<f64>::zeros(m, n);
63    for j in 0..n {
64        let mut norm_sq = 0.0;
65        for i in 0..m {
66            let v = mixing[(i, j)];
67            norm_sq += v * v;
68        }
69        let denom = norm_sq.sqrt().max(f64::EPSILON);
70        for i in 0..m {
71            au[(i, j)] = mixing[(i, j)] / denom;
72        }
73    }
74    let mut s = Mat::<f64>::zeros(n, n);
75    for i in 0..n {
76        for j in 0..n {
77            let mut dot = 0.0;
78            for r in 0..m {
79                dot += au[(r, i)] * au[(r, j)];
80            }
81            s[(i, j)] = dot;
82        }
83    }
84    hotspot_from_similarity(s.as_ref())
85}
86
87#[cfg(test)]
88mod tests {
89    use super::*;
90    use faer::Mat;
91
92    #[test]
93    fn identity_similarity_has_unit_sifs() {
94        let s = Mat::<f64>::from_fn(3, 3, |i, j| if i == j { 1.0 } else { 0.0 });
95        let h = hotspot_from_similarity(s.as_ref()).expect("hotspot");
96        for sif in h.sifs() {
97            assert!((sif - 1.0).abs() < 1e-9, "sif={sif}");
98        }
99    }
100
101    #[test]
102    fn collinear_columns_inflate_sifs() {
103        // Two nearly identical unit columns → high CS → large SIFs.
104        let a = Mat::<f64>::from_fn(4, 2, |i, j| {
105            if j == 0 {
106                if i == 0 {
107                    1.0
108                } else {
109                    0.05
110                }
111            } else if i == 0 {
112                0.98
113            } else if i == 1 {
114                0.2
115            } else {
116                0.05
117            }
118        });
119        let h = hotspot_from_mixing_matrix(a.as_ref()).expect("hotspot");
120        let sifs = h.sifs();
121        assert!(sifs[0] > 1.5, "sif0={}", sifs[0]);
122        assert!(sifs[1] > 1.5, "sif1={}", sifs[1]);
123        assert!(h.matrix[(0, 1)] > 0.5, "offdiag={}", h.matrix[(0, 1)]);
124    }
125
126    #[test]
127    fn singular_similarity_errors() {
128        let s = Mat::<f64>::from_fn(2, 2, |_i, _j| 1.0);
129        assert!(hotspot_from_similarity(s.as_ref()).is_err());
130    }
131}
flow_linalg/hotspot.rs

flow_linalg/
hotspot.rs
