qualia_core_db/solvers/learning/dimensionality/
som.rs

1//! Self-Organizing Map (CI-SKM ch 3) — a topology-preserving projection of a
2//! high-dimensional space onto a 2-D grid of neurons. Nearby inputs map to nearby
3//! grid cells, so it lays out a semantic space for the 10D→5D relevance router
4//! (complementing PCA: SOM preserves neighbourhood topology, not just variance).
5//! Kernel-class `AllPairs` (the best-matching-unit search). Deterministic per seed.
6
7use crate::solvers::learning::LearningError;
8
9/// A trained self-organizing map: a `grid_w × grid_h` lattice of `dim`-vectors.
10#[derive(Debug, Clone)]
11pub struct Som {
12    pub grid_w: usize,
13    pub grid_h: usize,
14    pub dim: usize,
15    weights: Vec<f64>, // grid_w*grid_h*dim, row-major over (y, x, d)
16}
17
18struct Lcg(u64);
19impl Lcg {
20    fn unit(&mut self) -> f64 {
21        self.0 = self
22            .0
23            .wrapping_mul(6364136223846793005)
24            .wrapping_add(1442695040888963407);
25        ((self.0 >> 11) as f64) / ((1u64 << 53) as f64)
26    }
27}
28
29impl Som {
30    fn idx(&self, x: usize, y: usize) -> usize {
31        (y * self.grid_w + x) * self.dim
32    }
33
34    /// Best-matching unit (grid coords) for an input vector — the nearest neuron.
35    pub fn bmu(&self, input: &[f64]) -> (usize, usize) {
36        let mut best = (0, 0);
37        let mut best_d = f64::INFINITY;
38        for y in 0..self.grid_h {
39            for x in 0..self.grid_w {
40                let w = &self.weights[self.idx(x, y)..self.idx(x, y) + self.dim];
41                let d: f64 = w.iter().zip(input).map(|(a, b)| (a - b) * (a - b)).sum();
42                if d < best_d {
43                    best_d = d;
44                    best = (x, y);
45                }
46            }
47        }
48        best
49    }
50
51    /// Train a `grid_w × grid_h` SOM on a row-major `n × dim` data matrix for
52    /// `epochs`, with initial learning rate `lr0` and neighbourhood radius `sigma0`
53    /// (both decaying exponentially). Fails closed on shape mismatch.
54    pub fn train(
55        data: &[f64],
56        n: usize,
57        dim: usize,
58        grid_w: usize,
59        grid_h: usize,
60        epochs: usize,
61        lr0: f64,
62        sigma0: f64,
63        seed: u64,
64    ) -> Result<Self, LearningError> {
65        if n == 0 || dim == 0 || grid_w == 0 || grid_h == 0 || data.len() != n * dim {
66            return Err(LearningError::InvalidDimension);
67        }
68        let mut rng = Lcg(seed ^ 0x9E3779B97F4A7C15);
69        // Initialise weights from the data range.
70        let mut lo = vec![f64::INFINITY; dim];
71        let mut hi = vec![f64::NEG_INFINITY; dim];
72        for i in 0..n {
73            for d in 0..dim {
74                let v = data[i * dim + d];
75                lo[d] = lo[d].min(v);
76                hi[d] = hi[d].max(v);
77            }
78        }
79        let mut weights = vec![0.0; grid_w * grid_h * dim];
80        for cell in 0..grid_w * grid_h {
81            for d in 0..dim {
82                weights[cell * dim + d] = lo[d] + rng.unit() * (hi[d] - lo[d]).max(1e-9);
83            }
84        }
85        let mut som = Som {
86            grid_w,
87            grid_h,
88            dim,
89            weights,
90        };
91
92        let total = epochs.max(1) as f64;
93        for epoch in 0..epochs.max(1) {
94            let frac = epoch as f64 / total;
95            let lr = lr0 * (-frac * 3.0).exp();
96            let sigma = (sigma0 * (-frac * 3.0).exp()).max(0.5);
97            let two_sigma2 = 2.0 * sigma * sigma;
98            for i in 0..n {
99                let input = &data[i * dim..(i + 1) * dim];
100                let (bx, by) = som.bmu(input);
101                for y in 0..grid_h {
102                    for x in 0..grid_w {
103                        let gd2 =
104                            ((x as f64 - bx as f64).powi(2)) + ((y as f64 - by as f64).powi(2));
105                        let h = (-gd2 / two_sigma2).exp();
106                        if h < 1e-4 {
107                            continue;
108                        }
109                        let base = som.idx(x, y);
110                        for d in 0..dim {
111                            som.weights[base + d] += lr * h * (input[d] - som.weights[base + d]);
112                        }
113                    }
114                }
115            }
116        }
117        Ok(som)
118    }
119}
120
121#[cfg(test)]
122mod tests {
123    use super::*;
124
125    #[test]
126    fn separated_clusters_map_to_distinct_regions() {
127        // Two clusters in 2-D, far apart; their BMUs should land in different grid
128        // regions (grid distance > 0), and same-cluster points should be close.
129        let mut data = Vec::new();
130        for d in 0..10 {
131            data.push((d % 3) as f64 * 0.1);
132            data.push((d % 3) as f64 * 0.1); // cluster near (0,0)
133        }
134        for d in 0..10 {
135            data.push(10.0 + (d % 3) as f64 * 0.1);
136            data.push(10.0 + (d % 3) as f64 * 0.1); // cluster near (10,10)
137        }
138        let som = Som::train(&data, 20, 2, 4, 4, 200, 0.5, 2.0, 1).unwrap();
139        let a = som.bmu(&[0.0, 0.0]);
140        let b = som.bmu(&[10.0, 10.0]);
141        let grid_dist =
142            ((a.0 as f64 - b.0 as f64).powi(2) + (a.1 as f64 - b.1 as f64).powi(2)).sqrt();
143        assert!(
144            grid_dist > 1.0,
145            "clusters should map apart: a={a:?} b={b:?}"
146        );
147        // A point near cluster A maps to A's BMU (or adjacent).
148        let a2 = som.bmu(&[0.2, 0.2]);
149        let near = ((a.0 as f64 - a2.0 as f64).powi(2) + (a.1 as f64 - a2.1 as f64).powi(2)).sqrt();
150        assert!(near <= 1.5, "same-cluster points map near: {a:?} vs {a2:?}");
151    }
152
153    #[test]
154    fn preserves_1d_order() {
155        // Inputs increasing along a line → BMUs should be (weakly) monotone on the grid.
156        let data: Vec<f64> = (0..10).map(|i| i as f64).collect();
157        let som = Som::train(&data, 10, 1, 10, 1, 300, 0.5, 3.0, 7).unwrap();
158        let first = som.bmu(&[0.0]).0 as i64;
159        let last = som.bmu(&[9.0]).0 as i64;
160        assert_ne!(first, last, "endpoints should occupy different cells");
161    }
162
163    #[test]
164    fn guards() {
165        assert_eq!(
166            Som::train(&[1.0], 1, 2, 3, 3, 10, 0.5, 1.0, 0).unwrap_err(),
167            LearningError::InvalidDimension
168        );
169    }
170}
qualia_core_db/solvers/learning/dimensionality/som.rs

qualia_core_db/solvers/learning/dimensionality/
som.rs
