flow_pacmap/
pca.rs

1//! PCA initialisation for PaCMAP.
2//!
3//! A two-component specialization of [`flow_dimensional_reduction::Pca`], which
4//! uses the covariance method: O(n·d²) to build a d×d matrix plus O(d³) to
5//! decompose it, rather than an O(n²)-scale decomposition of the data matrix.
6
7use crate::error::PaCMAPError;
8use flow_dimensional_reduction::Pca;
9
10/// Project `n × d` row-major data onto its top 2 principal components.
11///
12/// Returns `n` (PC1, PC2) score pairs.
13///
14/// # Errors
15/// Returns [`PaCMAPError::Pca`] if the decomposition fails or the inputs are
16/// inconsistent.
17pub fn pca_init(data: &[f32], n: usize, d: usize) -> Result<Vec<[f32; 2]>, PaCMAPError> {
18    debug_assert_eq!(data.len(), n * d);
19
20    let pca = Pca::new(2).fit(data, n, d)?;
21    let flat = pca.transform(data, n, d)?;
22    let k = pca.n_components();
23
24    // `k` is 1 when d == 1; pad the second axis with zeros to keep the
25    // [f32; 2] contract that callers rely on.
26    Ok(flat
27        .chunks(k)
28        .map(|c| [c[0], if k >= 2 { c[1] } else { 0.0 }])
29        .collect())
30}
31
32#[cfg(test)]
33mod tests {
34    use super::*;
35
36    #[test]
37    fn pca_separates_axis_aligned_clusters() {
38        // Two clusters separated along dim 0; PCA should put that on PC1
39        let mut data: Vec<f32> = Vec::new();
40        for _ in 0..50 {
41            data.extend_from_slice(&[0.0_f32, 0.0]);
42        }
43        for _ in 0..50 {
44            data.extend_from_slice(&[10.0_f32, 0.0]);
45        }
46        let emb = pca_init(&data, 100, 2).unwrap();
47        // First 50 points should have similar x; second 50 should differ
48        let mean_left_x = emb[..50].iter().map(|p| p[0]).sum::<f32>() / 50.0;
49        let mean_right_x = emb[50..].iter().map(|p| p[0]).sum::<f32>() / 50.0;
50        assert!(
51            (mean_left_x - mean_right_x).abs() > 1.0,
52            "PCA should separate the two clusters along PC1"
53        );
54    }
55}
flow_pacmap/pca.rs

flow_pacmap/
pca.rs
