Added Clustering, PCA, Sparse utilities

Author: Ian

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "single_algebra"
-version = "0.1.2-alpha.7"
+version = "0.2.0-alpha.0"
 edition = "2021"
 license-file = "LICENSE.md"
 description = "A linear algebra convenience library for the single-rust library. Can be used externally as well."
@@ -22,28 +22,36 @@ lapack = ["dep:nalgebra-lapack", "nalgebra-lapack/openblas", "dep:nalgebra"]
 nalgebra = ["dep:nalgebra"]
 faer = ["dep:faer", "dep:faer-ext"]
 simba = ["dep:simba"]
+clustering = ["network", "local_moving", "dep:kiddo"]
+network = []
+local_moving = ["network", "dep:ahash"]
 
 
 [dependencies]
 anyhow = "1.0.95"
 
 env_logger = "0.11.6"
-faer = {version = "0.20.1", optional = true}
-faer-ext = {version = "0.4.1", features = ["ndarray"], optional = true}
+faer = { version = "0.20.1", optional = true }
+faer-ext = { version = "0.4.1", features = ["ndarray"], optional = true }
 log = "0.4.22"
-nalgebra = {version = "0.33", features = ["serde-serialize"], optional = true}
-nalgebra-lapack = {version = "0.25.0", optional = true, default-features = false}
+nalgebra = { version = "0.33", features = ["serde-serialize"], optional = true }
+nalgebra-lapack = { version = "0.25.0", optional = true, default-features = false }
 nalgebra-sparse = "0.10"
-ndarray = {version = "0.16", features = ["rayon"]}
-nshare = {version = "0.10.0", features = ["ndarray", "nalgebra"]}
+ndarray = { version = "0.16", features = ["rayon"] }
+nshare = { version = "0.10.0", features = ["ndarray", "nalgebra"] }
 num-traits = "0.2.19"
 rayon = "1.10.0"
-simba = {version = "0.9.0", optional = true}
-smartcore = {version = "0.4", features=["ndarray-bindings"], optional = true}
+simba = { version = "0.9.0", optional = true }
+smartcore = { version = "0.4", features = ["ndarray-bindings"], optional = true }
 single-svdlib = "0.1.0"
+parking_lot = "0.12.3"
+petgraph = { version = "0.7.1", features = ["rayon"] }
+rand = "0.9.0"
+rand_chacha = "0.9.0"
+kiddo = { version = "5.0.3", optional = true }
+ahash = { version = "0.8.11", optional = true, features = ["compile-time-rng"] }
 
 
 [dev-dependencies]
 criterion = { version = "0.5.1", features = ["html_reports"] }
-rand = "0.8.5"
 approx = "0.5.1"

Cargo.toml

@@ -0,0 +1,41 @@
+# single-algebra 🧮
+
+The companion algebra library for single-rust, providing powerful matrix operations and machine learning utilities.
+
+## Features 🚀
+
+- Efficient operations on sparse and dense matrices
+- Dimensionality reduction techniques
+- Clustering algorithms including Louvain community detection
+- More features planned!
+
+## Matrix Operations 📊
+
+- SVD decomposition with parallel and LAPACK implementations
+- Matrix convenience functions for statistical operations
+- Support for both CSR and CSC sparse matrix formats
+
+## Clustering 🔍
+
+- Louvain community detection
+- Similarity network construction
+- K-nearest neighbors graph building
+- Local moving algorithm for community refinement
+
+## Dimensionality Reduction ⬇️
+
+- Incremental PCA implementation
+- Support for sparse matrices in dimensionality reduction
+
+## Acknowledgments 🙏
+
+The Louvain clustering implementation was adapted from [louvain-rs](https://github.com/graphext/louvain-rs/tree/master) written by Juan Morales ([email protected]). The original implementation has been modified to better suit the needs of single-algebra.
+
+## Installation
+
+Add this to your `Cargo.toml`:
+
+```toml
+[dependencies]
+single-algebra = "0.2.0-alpha.0"
+```

README.md

@@ -1,9 +1,9 @@
 use criterion::{BenchmarkId, Criterion, BenchmarkGroup, criterion_group, criterion_main};
 use nalgebra_sparse::{CooMatrix, CscMatrix};
-use rand::{distributions::Uniform, SeedableRng, rngs::StdRng};
+use rand::{SeedableRng, rngs::StdRng};
 use std::time::Duration;
 use criterion::measurement::Measurement;
-use rand::distributions::Distribution;
+use rand::distr::{Distribution, Uniform};
 use rayon::ThreadPool;
 use single_algebra::sparse::{MatrixNonZero, MatrixSum};
 
@@ -43,9 +43,9 @@ fn create_csc_matrix(
     let mut rng = StdRng::seed_from_u64(seed);
     let mut coo = CooMatrix::new(rows, cols);
     let total_elements = (rows * cols) as f64 * density;
-    let value_dist = Uniform::from(0.0..1.0);
-    let row_dist = Uniform::from(0..rows);
-    let col_dist = Uniform::from(0..cols);
+    let value_dist = Uniform::try_from(0.0..1.0).unwrap();
+    let row_dist = Uniform::try_from(0..rows).unwrap();
+    let col_dist = Uniform::try_from(0..cols).unwrap();
 
     for _ in 0..total_elements as usize {
         let row = row_dist.sample(&mut rng);

benches/csc_matrix_benchmark.rs

@@ -1,10 +1,10 @@
 use criterion::measurement::Measurement;
 use criterion::{criterion_group, criterion_main, BenchmarkGroup, BenchmarkId, Criterion};
 use nalgebra_sparse::{CooMatrix, CsrMatrix};
-use rand::distributions::Distribution;
-use rand::{distributions::Uniform, rngs::StdRng, SeedableRng};
+use rand::{rngs::StdRng, SeedableRng};
 use single_algebra::sparse::{MatrixNonZero, MatrixSum};
 use std::time::Duration;
+use rand::distr::{Distribution, Uniform};
 
 #[derive(Clone)]
 pub struct SparseMatrixConfig {
@@ -37,9 +37,9 @@ fn create_test_matrix(rows: usize, cols: usize, density: f64, seed: u64) -> CsrM
     let mut rng = StdRng::seed_from_u64(seed);
     let mut coo = CooMatrix::new(rows, cols);
     let total_elements = (rows * cols) as f64 * density;
-    let value_dist = Uniform::from(0.0..1.0);
-    let row_dist = Uniform::from(0..rows);
-    let col_dist = Uniform::from(0..cols);
+    let value_dist = Uniform::try_from(0.0..1.0).unwrap();
+    let row_dist = Uniform::try_from(0..rows).unwrap();
+    let col_dist = Uniform::try_from(0..cols).unwrap();
 
     for _ in 0..total_elements as usize {
         let row = row_dist.sample(&mut rng);

benches/csr_matrix_benchmark.rs

@@ -0,0 +1,172 @@
+// See https://www.nature.com/articles/s41598-019-41695-z#Sec2 and for the original publication: https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008 - https://en.wikipedia.org/wiki/Louvain_method & https://github.com/graphext/louvain-rs/tree/master
+// Copyright 2018 Juan Morales ([email protected])
+// Repository: https://github.com/graphext/louvain-rs/tree/master
+// Licensed under the MIT License.
+use std::collections::HashSet;
+use num_traits::{Float, FromPrimitive, ToPrimitive};
+use std::iter::Sum;
+use std::ops::MulAssign;
+use rand_chacha::ChaCha20Rng;
+use rand_chacha::rand_core::SeedableRng;
+use crate::local_moving::standard::StandardLocalMoving;
+use crate::network::{Graph, Network};
+use crate::network::clustering::{NetworkGrouping, VectorGrouping};
+
+pub const DEF_RES: f64 = 1.0;
+
+pub struct Louvain<T>
+where
+    T: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign, {
+    rng: ChaCha20Rng,
+    local_moving: StandardLocalMoving<T>
+}
+
+impl<T> Louvain<T>
+where
+    T: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign, {
+    pub fn new(resolution: T, seed: Option<u64>) -> Self {
+        let seed = seed.unwrap_or_default();
+
+        Louvain {
+            rng: ChaCha20Rng::seed_from_u64(seed),
+            local_moving: StandardLocalMoving::new(resolution)
+        }
+    }
+
+    pub fn iterate_one_level<N, E>(
+        &mut self,
+        network: &Network<N, E>,
+        clustering: &mut VectorGrouping
+    ) -> bool
+    where
+    N: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign,
+    E: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign, {
+        self.local_moving.iterate(network, clustering, &mut self.rng)
+    }
+
+    pub fn iterate<N, E>(
+        &mut self,
+        network: &Network<N, E>,
+        clustering: &mut VectorGrouping
+    ) -> bool
+    where
+        N: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign,
+        E: Float + FromPrimitive + ToPrimitive + Send + Sync + Sum + MulAssign {
+        let mut update = self.local_moving.iterate(network, clustering, &mut self.rng);
+
+        if clustering.group_count() == network.nodes() {
+            return update;
+        }
+
+        let reduced_network = network.create_reduced_network(clustering);
+        let mut reduced_clustering = VectorGrouping::create_isolated(reduced_network.nodes());
+        update |= self.iterate(&reduced_network, &mut reduced_clustering);
+        clustering.merge(&reduced_clustering);
+        update
+    }
+
+    pub fn build_network<I>(
+        n_nodes: usize,
+        n_edges: usize,
+        adjacency: I
+    ) -> Network<f64, f64>
+    where
+        I: Iterator<Item = (u32, u32)> {
+        let mut graph = Graph::with_capacity(n_nodes, n_edges);
+        let mut node_indices = Vec::with_capacity(n_nodes);
+
+        for _ in 0..n_nodes {
+            node_indices.push(graph.add_node(1.0));
+        }
+
+        let mut seen = vec![HashSet::<u32>::new(); n_nodes];
+        let mut node_weights = vec![0.0; n_nodes];
+
+        for (i, j) in adjacency {
+            let (i, j) = if i < j { (i, j) } else { (j, i) };
+            let i_ = i as usize;
+            let j_ = j as usize;
+
+            if seen[i_].insert(j) {
+                graph.add_edge(
+                    node_indices[i_],
+                    node_indices[j_],
+                    1.0
+                );
+                node_weights[j_] += 1.0;
+                node_weights[i_] += 1.0;
+            }
+        }
+
+        for &i in &node_indices {
+            *graph.node_weight_mut(i).unwrap() = node_weights[i.index()];
+        }
+
+        Network::new_from_graph(graph)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use petgraph::graph::NodeIndex;
+
+    fn create_test_network() -> Network<f64, f64> {
+        let mut graph = Graph::new_undirected();
+
+        // Add 5 nodes
+        for _ in 0..5 {
+            graph.add_node(1.0);
+        }
+
+        // Add edges to create two communities
+        graph.add_edge(NodeIndex::new(0), NodeIndex::new(1), 1.0);
+        graph.add_edge(NodeIndex::new(1), NodeIndex::new(2), 1.0);
+        graph.add_edge(NodeIndex::new(0), NodeIndex::new(2), 1.0);
+        graph.add_edge(NodeIndex::new(3), NodeIndex::new(4), 1.0);
+
+        Network::new_from_graph(graph)
+    }
+
+    #[test]
+    fn test_louvain_clustering() {
+        let network = create_test_network();
+        let mut clustering = VectorGrouping::create_isolated(network.nodes());
+        let mut louvain: Louvain<f64> = Louvain::new(DEF_RES.into(), Some(42));
+
+        assert!(louvain.iterate(&network, &mut clustering));
+
+        // Should identify two communities
+        assert!(clustering.group_count() == 2);
+
+        // Nodes 0,1,2 should be in same cluster
+        let cluster1 = clustering.get_group(0);
+        assert_eq!(clustering.get_group(1), cluster1);
+        assert_eq!(clustering.get_group(2), cluster1);
+
+        // Nodes 3,4 should be in different cluster
+        let cluster2 = clustering.get_group(3);
+        assert_eq!(clustering.get_group(4), cluster2);
+        assert_ne!(cluster1, cluster2);
+    }
+
+    #[test]
+    fn test_build_network() {
+        let edges = vec![(0, 1), (1, 2), (2, 0), (3, 4)];
+        let network = Louvain::<f64>::build_network(5, edges.len(), edges.into_iter());
+
+        assert_eq!(network.nodes(), 5);
+        assert_eq!(network.graph.edge_count(), 4);
+
+        // Check node weights (should equal degree)
+        for i in 0..5 {
+            let weight = network.weight(i);
+            let expected = match i {
+                0..=2 => 2.0,  // Nodes in triangle
+                3..=4 => 1.0,  // Nodes in single edge
+                _ => unreachable!(),
+            };
+            assert_eq!(weight, expected);
+        }
+    }
+}

src/pca/incremental/mod.rs src/clustering/leiden/mod.rssrc/pca/incremental/mod.rs renamed to src/clustering/leiden/mod.rs

@@ -0,0 +1,7 @@
+pub(crate) mod leiden;
+pub(crate) mod louvain;
+pub use louvain::Louvain;
+pub(crate) mod similarity_network;
+pub use similarity_network::build_knn_network_combined_matrix;
+pub use similarity_network::build_knn_network_separate_matrix;
+pub use similarity_network::create_similarity_network;