anofox_ml_manifold/
lib.rs

1//! Manifold learning algorithms.
2//!
3//! Provides:
4//! - **Classical MDS** — eigendecomposition variant of `sklearn.manifold.MDS`.
5//! - **Isomap** — `sklearn.manifold.Isomap`: k-NN graph + geodesic shortest
6//!   paths (Floyd-Warshall) + classical MDS on the geodesic distances.
7//!
8//! - **LocallyLinearEmbedding** — `sklearn.manifold.LocallyLinearEmbedding`:
9//!   local reconstruction weights + bottom-k eigenvectors of `(I − W)ᵀ(I − W)`.
10//!
11//! - **t-SNE** — `sklearn.manifold.TSNE`: vanilla O(n²) gradient descent
12//!   with perplexity-calibrated Gaussian affinities and student-t low-dim
13//!   affinities. Suitable for n ≲ 1000.
14
15pub mod isomap;
16pub mod lle;
17pub mod tsne;
18
19pub use isomap::{FittedIsomap, Isomap};
20pub use lle::{FittedLocallyLinearEmbedding, LocallyLinearEmbedding};
21pub use tsne::{FittedTSne, TSne, TSneMethod};
22
23use anofox_ml_core::{FitUnsupervised, Result, RustMlError};
24use faer::linalg::solvers::SelfAdjointEigen;
25use faer::{Mat, Side};
26use ndarray::{Array1, Array2};
27
28#[derive(Debug, Clone)]
29pub struct ClassicalMds {
30    pub n_components: usize,
31}
32
33impl ClassicalMds {
34    pub fn new(n_components: usize) -> Self {
35        Self { n_components }
36    }
37}
38
39#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
40pub struct FittedClassicalMds {
41    pub embedding: Array2<f64>,
42    pub eigenvalues: Array1<f64>,
43}
44
45fn pairwise_dist(x: &Array2<f64>) -> Array2<f64> {
46    let n = x.nrows();
47    let mut d = Array2::<f64>::zeros((n, n));
48    for i in 0..n {
49        let xi = x.row(i).to_owned();
50        for j in (i + 1)..n {
51            let xj = x.row(j).to_owned();
52            let mut sd = 0.0;
53            for k in 0..x.ncols() {
54                sd += (xi[k] - xj[k]).powi(2);
55            }
56            let v = sd.sqrt();
57            d[[i, j]] = v;
58            d[[j, i]] = v;
59        }
60    }
61    d
62}
63
64impl FitUnsupervised<f64> for ClassicalMds {
65    type Fitted = FittedClassicalMds;
66
67    fn fit(&self, x: &Array2<f64>) -> Result<Self::Fitted> {
68        let n = x.nrows();
69        let k = self.n_components.min(n);
70        if k == 0 {
71            return Err(RustMlError::InvalidParameter("n_components >= 1".into()));
72        }
73        let d = pairwise_dist(x);
74
75        // Double-centered squared distance: B = -1/2 * (I - 1/n J) * D² * (I - 1/n J)
76        let mut d2 = Array2::<f64>::zeros((n, n));
77        for i in 0..n {
78            for j in 0..n {
79                d2[[i, j]] = d[[i, j]] * d[[i, j]];
80            }
81        }
82        let mut row_mean = vec![0.0_f64; n];
83        let mut col_mean = vec![0.0_f64; n];
84        let mut global = 0.0_f64;
85        for i in 0..n {
86            for j in 0..n {
87                row_mean[i] += d2[[i, j]];
88                col_mean[j] += d2[[i, j]];
89                global += d2[[i, j]];
90            }
91        }
92        let n_f = n as f64;
93        for i in 0..n {
94            row_mean[i] /= n_f;
95            col_mean[i] /= n_f;
96        }
97        let global = global / (n_f * n_f);
98
99        let mut b = Array2::<f64>::zeros((n, n));
100        for i in 0..n {
101            for j in 0..n {
102                b[[i, j]] = -0.5 * (d2[[i, j]] - row_mean[i] - col_mean[j] + global);
103            }
104        }
105
106        // Symmetric eigendecomposition; eigenvalues come back ascending.
107        let m = Mat::<f64>::from_fn(n, n, |i, j| b[[i, j]]);
108        let eig = SelfAdjointEigen::new(m.as_ref(), Side::Lower)
109            .map_err(|e| RustMlError::InvalidParameter(format!("eigen failed: {e:?}")))?;
110        let s = eig.S();
111        let v = eig.U();
112
113        let mut embedding = Array2::<f64>::zeros((n, k));
114        let mut eigenvalues = Array1::<f64>::zeros(k);
115        for c in 0..k {
116            let src = n - 1 - c;
117            let lam = s.column_vector()[src];
118            eigenvalues[c] = lam;
119            let scale = lam.max(0.0).sqrt();
120            for i in 0..n {
121                embedding[[i, c]] = v[(i, src)] * scale;
122            }
123        }
124        Ok(FittedClassicalMds {
125            embedding,
126            eigenvalues,
127        })
128    }
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134    use ndarray::array;
135
136    #[test]
137    fn test_mds_recovers_planar_distances() {
138        // Points in 2D — MDS to 2 components should preserve pairwise distances.
139        let x = array![
140            [0.0_f64, 0.0],
141            [1.0, 0.0],
142            [0.0, 1.0],
143            [2.0, 2.0],
144            [-1.0, 3.0],
145        ];
146        let mds = ClassicalMds::new(2);
147        let fitted = mds.fit(&x).unwrap();
148        // Pairwise distances in the embedding should match the originals.
149        let orig = pairwise_dist(&x);
150        let emb = pairwise_dist(&fitted.embedding);
151        for i in 0..5 {
152            for j in 0..5 {
153                assert!(
154                    (orig[[i, j]] - emb[[i, j]]).abs() < 1e-6,
155                    "dist[{i},{j}]: orig {} vs emb {}",
156                    orig[[i, j]],
157                    emb[[i, j]]
158                );
159            }
160        }
161    }
162}
anofox_ml_manifold/lib.rs

anofox_ml_manifold/
lib.rs
