From ee83e62e94ec1bdb708fbfd8681d4b0844828bb5 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 21:21:24 +0100
Subject: [PATCH 01/22] Scaffolding

---
 Cargo.toml          |  1 +
 k_means/Cargo.toml  | 19 +++++++++++++++++++
 k_means/src/lib.rs  |  0
 k_means/src/main.rs |  3 +++
 4 files changed, 23 insertions(+)
 create mode 100644 k_means/Cargo.toml
 create mode 100644 k_means/src/lib.rs
 create mode 100644 k_means/src/main.rs

diff --git a/Cargo.toml b/Cargo.toml
index 54b5fe5..9335a04 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -2,4 +2,5 @@
 members = [
   "./",
   "linear_regression",
+  "k_means",
 ]
\ No newline at end of file
diff --git a/k_means/Cargo.toml b/k_means/Cargo.toml
new file mode 100644
index 0000000..6a4ac9e
--- /dev/null
+++ b/k_means/Cargo.toml
@@ -0,0 +1,19 @@
+[package]
+name = "k_means"
+version = "0.1.0"
+authors = ["LukeMathWalker"]
+edition = "2018"
+workspace = ".."
+
+[features]
+default = []
+openblas = ["ndarray-linalg/openblas"]
+intel-mkl = ["ndarray-linalg/intel-mkl"]
+netlib = ["ndarray-linalg/netlib"]
+
+[dependencies]
+ndarray = {version = "0.12", features = ["blas"]}
+ndarray-linalg = {version = "0.11.1", optional = true, default-features = false}
+ndarray-stats = {git = "https://github.com/rust-ndarray/ndarray-stats", branch = "master"}
+ndarray-rand = "0.9"
+rand = "0.6"
diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
new file mode 100644
index 0000000..e69de29
diff --git a/k_means/src/main.rs b/k_means/src/main.rs
new file mode 100644
index 0000000..9d82917
--- /dev/null
+++ b/k_means/src/main.rs
@@ -0,0 +1,3 @@
+pub fn main() {
+
+}
\ No newline at end of file

From 4e4692af00d86eda5f0e2f731a6062689977d7a9 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 22:08:00 +0100
Subject: [PATCH 02/22] Started the implementation. Blocked by ndarray-stats
 issue on data ownership for square l2 distance

---
 k_means/src/lib.rs  | 67 +++++++++++++++++++++++++++++++++++++++++++++
 k_means/src/main.rs |  4 +--
 2 files changed, 69 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index e69de29..ae7531d 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -0,0 +1,67 @@
+#![allow(non_snake_case)]
+use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix2, Ix1};
+use rand::distributions::{Distribution, Uniform};
+use ndarray_stats::DeviationExt;
+
+pub struct KMeans {
+    pub centroids: Option<Array2<f64>>,
+    n_clusters: u16,
+}
+
+impl KMeans {
+    pub fn new(n_clusters: u16) -> KMeans {
+        KMeans {
+            centroids: None,
+            n_clusters,
+        }
+    }
+
+    /// Given an input matrix `X`, with shape `(n_samples, n_features)`
+    /// `fit` determines a set of centroids based on the training data distribution.
+    ///
+    /// `self` is modified in place, nothing is returned.
+    pub fn fit<A, B>(&mut self, X: &ArrayBase<A, Ix2>)
+    where
+        A: Data<Elem = f64>,
+    {
+        let centroids = KMeans::get_random_centroids(self.n_clusters, X);
+    }
+
+    fn get_random_centroids<A>(n_clusters: u16, X: &ArrayBase<A, Ix2>) -> Array2<f64>
+    where
+        A: Data<Elem = f64>,
+    {
+        let (n_samples, _) = X.dim();
+        let distribution = Uniform::from(0..n_samples);
+        let mut rng = rand::thread_rng();
+        let indices: Vec<_> = (0..n_clusters)
+            .map(|_| distribution.sample(&mut rng))
+            .collect();
+        X.select(Axis(0), &indices)
+    }
+
+    fn find_closest_centroid<A, B>(
+        centroids: &ArrayBase<A, Ix2>,
+        sample: &ArrayBase<B, Ix1>,
+    ) -> usize
+    where
+        A: Data<Elem = f64>,
+        B: Data<Elem = f64>,
+    {
+        let iterator = centroids.genrows().into_iter().enumerate();
+
+        let (_, first_centroid) = iterator.next().expect("No centroids - degenerate case!");
+        let mut closest_index = 0;
+        let mut minimum_distance = sample.sq_l2_dist(&first_centroid).unwrap();
+
+        for (index, centroid) in iterator {
+            let distance = sample.sq_l2_dist(&centroid).unwrap();
+            if distance < minimum_distance {
+                minimum_distance = distance;
+                closest_index = index;
+            }
+        }
+
+        closest_index
+    }
+}
diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index 9d82917..e842174 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -1,3 +1,3 @@
-pub fn main() {
+use k_means::KMeans;
 
-}
\ No newline at end of file
+pub fn main() {}

From 02154561c7b0afe7ab7bae807fb38d3ab7ecb0b2 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 22:26:20 +0100
Subject: [PATCH 03/22] Point to fork of ndarray-stats to make code compile

---
 k_means/Cargo.toml | 2 +-
 k_means/src/lib.rs | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/k_means/Cargo.toml b/k_means/Cargo.toml
index 6a4ac9e..ef092ad 100644
--- a/k_means/Cargo.toml
+++ b/k_means/Cargo.toml
@@ -14,6 +14,6 @@ netlib = ["ndarray-linalg/netlib"]
 [dependencies]
 ndarray = {version = "0.12", features = ["blas"]}
 ndarray-linalg = {version = "0.11.1", optional = true, default-features = false}
-ndarray-stats = {git = "https://github.com/rust-ndarray/ndarray-stats", branch = "master"}
+ndarray-stats = {git = "https://github.com/LukeMathWalker/ndarray-stats", branch = "generalized-deviation"}
 ndarray-rand = "0.9"
 rand = "0.6"
diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index ae7531d..c601a18 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -1,5 +1,5 @@
 #![allow(non_snake_case)]
-use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix2, Ix1};
+use ndarray::{Array2, ArrayBase, Axis, Data, Ix2, Ix1};
 use rand::distributions::{Distribution, Uniform};
 use ndarray_stats::DeviationExt;
 
@@ -48,7 +48,7 @@ impl KMeans {
         A: Data<Elem = f64>,
         B: Data<Elem = f64>,
     {
-        let iterator = centroids.genrows().into_iter().enumerate();
+        let mut iterator = centroids.genrows().into_iter().enumerate();
 
         let (_, first_centroid) = iterator.next().expect("No centroids - degenerate case!");
         let mut closest_index = 0;

From 11653902501bfa3bf1a23c519335cbdf4fe6c81b Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 22:26:56 +0100
Subject: [PATCH 04/22] Drop unneeded dependencies

---
 k_means/Cargo.toml | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/k_means/Cargo.toml b/k_means/Cargo.toml
index ef092ad..3a61d25 100644
--- a/k_means/Cargo.toml
+++ b/k_means/Cargo.toml
@@ -5,15 +5,8 @@ authors = ["LukeMathWalker"]
 edition = "2018"
 workspace = ".."
 
-[features]
-default = []
-openblas = ["ndarray-linalg/openblas"]
-intel-mkl = ["ndarray-linalg/intel-mkl"]
-netlib = ["ndarray-linalg/netlib"]
-
 [dependencies]
-ndarray = {version = "0.12", features = ["blas"]}
-ndarray-linalg = {version = "0.11.1", optional = true, default-features = false}
+ndarray = {version = "0.12"}
 ndarray-stats = {git = "https://github.com/LukeMathWalker/ndarray-stats", branch = "generalized-deviation"}
 ndarray-rand = "0.9"
 rand = "0.6"

From 1a45837b813735ac9ad94eb7b1b1b9cec3a9b346 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 22:31:59 +0100
Subject: [PATCH 05/22] Initialisation complete

---
 k_means/src/lib.rs | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index c601a18..76948c9 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -1,7 +1,7 @@
 #![allow(non_snake_case)]
-use ndarray::{Array2, ArrayBase, Axis, Data, Ix2, Ix1};
-use rand::distributions::{Distribution, Uniform};
+use ndarray::{Array2, ArrayBase, Axis, Data, Ix1, Ix2};
 use ndarray_stats::DeviationExt;
+use rand::distributions::{Distribution, Uniform};
 
 pub struct KMeans {
     pub centroids: Option<Array2<f64>>,
@@ -24,7 +24,11 @@ impl KMeans {
     where
         A: Data<Elem = f64>,
     {
+        // Initialisation
         let centroids = KMeans::get_random_centroids(self.n_clusters, X);
+        let cluster_memberships = X.map_axis(Axis(0), |sample| {
+            KMeans::find_closest_centroid(&centroids, &sample)
+        });
     }
 
     fn get_random_centroids<A>(n_clusters: u16, X: &ArrayBase<A, Ix2>) -> Array2<f64>

From 1def65318335263609019bd5a56c9a07b458677f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 6 Aug 2019 22:35:41 +0100
Subject: [PATCH 06/22] Add input check, clean iterator code

---
 k_means/src/lib.rs | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 76948c9..64e60b0 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -24,6 +24,12 @@ impl KMeans {
     where
         A: Data<Elem = f64>,
     {
+        let (n_samples, _) = X.dim();
+        assert!(
+            n_samples >= self.n_clusters as usize,
+            "We need more sample points than clusters!"
+        );
+
         // Initialisation
         let centroids = KMeans::get_random_centroids(self.n_clusters, X);
         let cluster_memberships = X.map_axis(Axis(0), |sample| {
@@ -52,13 +58,13 @@ impl KMeans {
         A: Data<Elem = f64>,
         B: Data<Elem = f64>,
     {
-        let mut iterator = centroids.genrows().into_iter().enumerate();
+        let mut iterator = centroids.genrows().into_iter();
 
-        let (_, first_centroid) = iterator.next().expect("No centroids - degenerate case!");
+        let first_centroid = iterator.next().expect("No centroids - degenerate case!");
         let mut closest_index = 0;
         let mut minimum_distance = sample.sq_l2_dist(&first_centroid).unwrap();
 
-        for (index, centroid) in iterator {
+        for (index, centroid) in iterator.enumerate() {
             let distance = sample.sq_l2_dist(&centroid).unwrap();
             if distance < minimum_distance {
                 minimum_distance = distance;

From 4a519793cccbc0ad4758b099655228776bb34f0c Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 14:48:21 +0100
Subject: [PATCH 07/22] Fix compilation errors

---
 k_means/src/lib.rs | 57 +++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 56 insertions(+), 1 deletion(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 64e60b0..b08c276 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -1,6 +1,7 @@
 #![allow(non_snake_case)]
-use ndarray::{Array2, ArrayBase, Axis, Data, Ix1, Ix2};
+use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2};
 use ndarray_stats::DeviationExt;
+use std::collections::HashMap;
 use rand::distributions::{Distribution, Uniform};
 
 pub struct KMeans {
@@ -35,6 +36,39 @@ impl KMeans {
         let cluster_memberships = X.map_axis(Axis(0), |sample| {
             KMeans::find_closest_centroid(&centroids, &sample)
         });
+
+        let centroids = KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
+    }
+
+    fn compute_centroids<A, B>(
+        X: &ArrayBase<A, Ix2>,
+        cluster_memberships: &ArrayBase<B, Ix1>,
+        n_clusters: u16
+    ) -> Array2<f64>
+    where
+        A: Data<Elem = f64>,
+        B: Data<Elem = usize>,
+    {
+        let (_, n_features) = X.dim();
+        let mut centroids: HashMap<usize, RollingMean> = HashMap::new();
+
+        let mut iterator = X.genrows().into_iter().zip(cluster_memberships.iter());
+        for (row, cluster_index) in iterator {
+            if let Some(rolling_mean) = centroids.get_mut(cluster_index) {
+                rolling_mean.accumulate(&row);
+            } else {
+                let new_centroid = RollingMean::new(row.to_owned());
+                centroids.insert(*cluster_index, new_centroid);
+            }
+        }
+
+        let mut new_centroids: Array2<f64> = Array2::zeros((n_clusters as usize, n_features));
+        for (cluster_index, centroid) in centroids.into_iter() {
+            let mut new_centroid = new_centroids.index_axis_mut(Axis(0), cluster_index);
+            new_centroid.assign(&centroid.current_mean);
+        }
+
+        new_centroids
     }
 
     fn get_random_centroids<A>(n_clusters: u16, X: &ArrayBase<A, Ix2>) -> Array2<f64>
@@ -75,3 +109,24 @@ impl KMeans {
         closest_index
     }
 }
+
+struct RollingMean {
+    pub current_mean: Array1<f64>,
+    n_samples: u64,
+}
+
+impl RollingMean {
+    pub fn new(first_sample: Array1<f64>) -> Self {
+        RollingMean { current_mean: first_sample, n_samples: 1 }
+    }
+
+    pub fn accumulate<A>(&mut self, new_sample: &ArrayBase<A, Ix1>)
+    where
+        A: Data<Elem = f64>,
+    {
+        let mut increment: Array1<f64> = &self.current_mean - new_sample;
+        increment.mapv_inplace(|x| x / (self.n_samples + 1) as f64);
+        self.current_mean -= &increment;
+        self.n_samples += 1;
+    }
+}

From d2ce0373f094ca71bb5dda205ac37e8c3ed41742 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 14:52:12 +0100
Subject: [PATCH 08/22] Update docs

---
 k_means/src/lib.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index b08c276..a2992f6 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -18,7 +18,7 @@ impl KMeans {
     }
 
     /// Given an input matrix `X`, with shape `(n_samples, n_features)`
-    /// `fit` determines a set of centroids based on the training data distribution.
+    /// `fit` determines `self.n_clusters` centroids based on the training data distribution.
     ///
     /// `self` is modified in place, nothing is returned.
     pub fn fit<A, B>(&mut self, X: &ArrayBase<A, Ix2>)
@@ -52,7 +52,7 @@ impl KMeans {
         let (_, n_features) = X.dim();
         let mut centroids: HashMap<usize, RollingMean> = HashMap::new();
 
-        let mut iterator = X.genrows().into_iter().zip(cluster_memberships.iter());
+        let iterator = X.genrows().into_iter().zip(cluster_memberships.iter());
         for (row, cluster_index) in iterator {
             if let Some(rolling_mean) = centroids.get_mut(cluster_index) {
                 rolling_mean.accumulate(&row);

From c3b80a09c18b9873cae1a88dfa1b214b428595da Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 14:59:03 +0100
Subject: [PATCH 09/22] Add training loop and naive convergence condition

---
 k_means/src/lib.rs | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index a2992f6..e3a7c49 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -31,13 +31,23 @@ impl KMeans {
             "We need more sample points than clusters!"
         );
 
+        let mut has_converged = false;
+
         // Initialisation
         let centroids = KMeans::get_random_centroids(self.n_clusters, X);
-        let cluster_memberships = X.map_axis(Axis(0), |sample| {
-            KMeans::find_closest_centroid(&centroids, &sample)
-        });
 
-        let centroids = KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
+        while !has_converged {
+            // Assignment step: associate each sample to the closest centroid
+            let cluster_memberships = X.map_axis(Axis(0), |sample| {
+                KMeans::find_closest_centroid(&centroids, &sample)
+            });
+
+            // Update step: calculate the mean of each cluster and use it as the new centroid
+            let new_centroids = KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
+
+            // Check convergence condition (very naive, we need an epsilon tolerance here)
+            has_converged = centroids == new_centroids;
+        }
     }
 
     fn compute_centroids<A, B>(

From 698ae8eafe16b8da079dc611d4da39f19ae90d4e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 15:07:16 +0100
Subject: [PATCH 10/22] Put together somethings that can be run

---
 k_means/src/lib.rs  |  4 ++--
 k_means/src/main.rs | 31 ++++++++++++++++++++++++++++++-
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index e3a7c49..736ada5 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -21,7 +21,7 @@ impl KMeans {
     /// `fit` determines `self.n_clusters` centroids based on the training data distribution.
     ///
     /// `self` is modified in place, nothing is returned.
-    pub fn fit<A, B>(&mut self, X: &ArrayBase<A, Ix2>)
+    pub fn fit<A>(&mut self, X: &ArrayBase<A, Ix2>)
     where
         A: Data<Elem = f64>,
     {
@@ -38,7 +38,7 @@ impl KMeans {
 
         while !has_converged {
             // Assignment step: associate each sample to the closest centroid
-            let cluster_memberships = X.map_axis(Axis(0), |sample| {
+            let cluster_memberships = X.map_axis(Axis(1), |sample| {
                 KMeans::find_closest_centroid(&centroids, &sample)
             });
 
diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index e842174..69da1cf 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -1,3 +1,32 @@
+#![allow(non_snake_case)]
+use ndarray::{Array, Array1, Array2, Axis};
+use ndarray_rand::RandomExt;
+use ndarray_stats::DeviationExt;
+use rand::distributions::{Normal, StandardNormal};
+
+// Import KMeans from other file ("lib.rs") in this example
 use k_means::KMeans;
 
-pub fn main() {}
+/// It returns a data distribution.
+///
+fn get_data(n_samples: usize, n_features: usize) -> Array2<f64> {
+    let shape = (n_samples, n_features);
+    let X: Array2<f64> = Array::random(shape, StandardNormal);
+    X
+}
+
+pub fn main() {
+    let n_samples = 5000;
+    let n_features = 3;
+    let n_clusters = 2;
+
+    let X = get_data(n_samples, n_features);
+
+    let mut k_means = KMeans::new(n_clusters);
+    k_means.fit(&X);
+
+    println!(
+        "The centroids are {:.3}",
+        k_means.centroids.unwrap()
+    );
+}

From 0ec1f3efa84919c1c975216517b2eb807637031d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 15:12:11 +0100
Subject: [PATCH 11/22] It works, more or less

---
 k_means/src/lib.rs | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 736ada5..25d1b12 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -32,9 +32,10 @@ impl KMeans {
         );
 
         let mut has_converged = false;
+        let tolerance = 1e-3;
 
         // Initialisation
-        let centroids = KMeans::get_random_centroids(self.n_clusters, X);
+        let mut centroids = KMeans::get_random_centroids(self.n_clusters, X);
 
         while !has_converged {
             // Assignment step: associate each sample to the closest centroid
@@ -46,8 +47,14 @@ impl KMeans {
             let new_centroids = KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
 
             // Check convergence condition (very naive, we need an epsilon tolerance here)
-            has_converged = centroids == new_centroids;
+            let distance = centroids.sq_l2_dist(&new_centroids).unwrap();
+            has_converged = distance < tolerance;
+            println!("Centroid distance: {:?}", distance);
+
+            centroids = new_centroids;
         }
+
+        self.centroids = Some(centroids);
     }
 
     fn compute_centroids<A, B>(

From 2da2b25c4dbd3e3dd0328768c66f3f565c3a15a4 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 17:37:02 +0100
Subject: [PATCH 12/22] Working!

---
 k_means/src/lib.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 25d1b12..e9dc56d 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -49,7 +49,6 @@ impl KMeans {
             // Check convergence condition (very naive, we need an epsilon tolerance here)
             let distance = centroids.sq_l2_dist(&new_centroids).unwrap();
             has_converged = distance < tolerance;
-            println!("Centroid distance: {:?}", distance);
 
             centroids = new_centroids;
         }
@@ -119,7 +118,8 @@ impl KMeans {
             let distance = sample.sq_l2_dist(&centroid).unwrap();
             if distance < minimum_distance {
                 minimum_distance = distance;
-                closest_index = index;
+                // We skipped the first centroid in the for loop
+                closest_index = index + 1;
             }
         }
 

From f4e0bd32d195617c93f80f3510472ff37593e7a1 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 17:56:27 +0100
Subject: [PATCH 13/22] Add comments

---
 k_means/src/main.rs | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index 69da1cf..48c0fbd 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -1,24 +1,27 @@
 #![allow(non_snake_case)]
-use ndarray::{Array, Array1, Array2, Axis};
+use ndarray::{Array, Array2, Axis, stack};
 use ndarray_rand::RandomExt;
-use ndarray_stats::DeviationExt;
-use rand::distributions::{Normal, StandardNormal};
+use rand::distributions::Normal;
 
 // Import KMeans from other file ("lib.rs") in this example
 use k_means::KMeans;
 
 /// It returns a data distribution.
 ///
+/// The data is clearly centered around two distinct points,
+/// to quickly spot if something is wrong with the KMeans algorithm
+/// looking at the output.
 fn get_data(n_samples: usize, n_features: usize) -> Array2<f64> {
-    let shape = (n_samples, n_features);
-    let X: Array2<f64> = Array::random(shape, StandardNormal);
-    X
+    let shape = (n_samples / 2, n_features);
+    let X1: Array2<f64> = Array::random(shape, Normal::new(1000., 0.1));
+    let X2: Array2<f64> = Array::random(shape, Normal::new(-1000., 0.1));
+    stack(Axis(0), &[X1.view(), X2.view()]).unwrap().to_owned()
 }
 
 pub fn main() {
-    let n_samples = 5000;
+    let n_samples = 5000000;
     let n_features = 3;
-    let n_clusters = 2;
+    let n_clusters = 8;
 
     let X = get_data(n_samples, n_features);
 

From dca79a01c6adfafdeff01155cfbf14be7ad7b4bc Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 17:59:55 +0100
Subject: [PATCH 14/22] Make it fast for the CI run

---
 k_means/src/main.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index 48c0fbd..4b987bc 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -19,9 +19,9 @@ fn get_data(n_samples: usize, n_features: usize) -> Array2<f64> {
 }
 
 pub fn main() {
-    let n_samples = 5000000;
+    let n_samples = 50000;
     let n_features = 3;
-    let n_clusters = 8;
+    let n_clusters = 2;
 
     let X = get_data(n_samples, n_features);
 

From 938522c96e8cb350ef534cb028a09a43af9533f3 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:00:19 +0100
Subject: [PATCH 15/22] Fmt

---
 k_means/src/lib.rs  | 12 ++++++++----
 k_means/src/main.rs |  7 ++-----
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index e9dc56d..b8e520f 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -1,8 +1,8 @@
 #![allow(non_snake_case)]
 use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2};
 use ndarray_stats::DeviationExt;
-use std::collections::HashMap;
 use rand::distributions::{Distribution, Uniform};
+use std::collections::HashMap;
 
 pub struct KMeans {
     pub centroids: Option<Array2<f64>>,
@@ -44,7 +44,8 @@ impl KMeans {
             });
 
             // Update step: calculate the mean of each cluster and use it as the new centroid
-            let new_centroids = KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
+            let new_centroids =
+                KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
 
             // Check convergence condition (very naive, we need an epsilon tolerance here)
             let distance = centroids.sq_l2_dist(&new_centroids).unwrap();
@@ -59,7 +60,7 @@ impl KMeans {
     fn compute_centroids<A, B>(
         X: &ArrayBase<A, Ix2>,
         cluster_memberships: &ArrayBase<B, Ix1>,
-        n_clusters: u16
+        n_clusters: u16,
     ) -> Array2<f64>
     where
         A: Data<Elem = f64>,
@@ -134,7 +135,10 @@ struct RollingMean {
 
 impl RollingMean {
     pub fn new(first_sample: Array1<f64>) -> Self {
-        RollingMean { current_mean: first_sample, n_samples: 1 }
+        RollingMean {
+            current_mean: first_sample,
+            n_samples: 1,
+        }
     }
 
     pub fn accumulate<A>(&mut self, new_sample: &ArrayBase<A, Ix1>)
diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index 4b987bc..500e19c 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -1,5 +1,5 @@
 #![allow(non_snake_case)]
-use ndarray::{Array, Array2, Axis, stack};
+use ndarray::{stack, Array, Array2, Axis};
 use ndarray_rand::RandomExt;
 use rand::distributions::Normal;
 
@@ -28,8 +28,5 @@ pub fn main() {
     let mut k_means = KMeans::new(n_clusters);
     k_means.fit(&X);
 
-    println!(
-        "The centroids are {:.3}",
-        k_means.centroids.unwrap()
-    );
+    println!("The centroids are {:.3}", k_means.centroids.unwrap());
 }

From bd9db67323449a43ed4976dc464f8cbaa1254e69 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:07:44 +0100
Subject: [PATCH 16/22] Adding comments

---
 k_means/src/lib.rs | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index b8e520f..63f0f8d 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -4,8 +4,28 @@ use ndarray_stats::DeviationExt;
 use rand::distributions::{Distribution, Uniform};
 use std::collections::HashMap;
 
+/// K-means clustering aims to partition a set of observations
+/// into `self.n_clusters` clusters, where each observation belongs
+/// to the cluster with the nearest mean.
+///
+/// The mean of the points within a cluster is called *centroid*.
+///
+/// Given the set of `centroids`, you can assign an observation to a cluster
+/// choosing the nearest centroid.
+///
+/// Details on the algorithm can be found [here](https://en.wikipedia.org/wiki/K-means_clustering).
+///
+/// We are implementing the _standard algorithm_.
 pub struct KMeans {
+    /// Our set of centroids.
+    ///
+    /// Before `fit` is called, it's set to `None`.
+    ///
+    /// Once `fit` is called, we will have our set of centroids: the `centroids` matrix
+    /// has shape `(n_clusters, n_features)`.
     pub centroids: Option<Array2<f64>>,
+    /// The number of clusters we are trying to subdivide our observations into.
+    /// It's set before-hand.
     n_clusters: u16,
 }
 

From c6b2e22a1e08678bee84aa84b9f685a0104c69ad Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:20:02 +0100
Subject: [PATCH 17/22] Add comments

---
 k_means/src/lib.rs | 25 ++++++++++++++++++-------
 1 file changed, 18 insertions(+), 7 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 63f0f8d..33c7f9f 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -30,6 +30,8 @@ pub struct KMeans {
 }
 
 impl KMeans {
+    /// The number of clusters we are looking for has to be chosen before-hand,
+    /// hence we take it as a constructor parameter of our `KMeans` struct.
     pub fn new(n_clusters: u16) -> KMeans {
         KMeans {
             centroids: None,
@@ -40,7 +42,7 @@ impl KMeans {
     /// Given an input matrix `X`, with shape `(n_samples, n_features)`
     /// `fit` determines `self.n_clusters` centroids based on the training data distribution.
     ///
-    /// `self` is modified in place, nothing is returned.
+    /// `self` is modified in place (`self.centroids` is mutated), nothing is returned.
     pub fn fit<A>(&mut self, X: &ArrayBase<A, Ix2>)
     where
         A: Data<Elem = f64>,
@@ -51,29 +53,38 @@ impl KMeans {
             "We need more sample points than clusters!"
         );
 
-        let mut has_converged = false;
         let tolerance = 1e-3;
 
-        // Initialisation
+        // Initialisation: we use the Forgy method - we pick `self.n_clusters` random
+        // observations and we use them as our first set of centroids.
         let mut centroids = KMeans::get_random_centroids(self.n_clusters, X);
 
-        while !has_converged {
+        // Keep repeating the assignment-update steps until we have convergence
+        loop {
             // Assignment step: associate each sample to the closest centroid
             let cluster_memberships = X.map_axis(Axis(1), |sample| {
                 KMeans::find_closest_centroid(&centroids, &sample)
             });
 
-            // Update step: calculate the mean of each cluster and use it as the new centroid
+            // Update step: using the newly computed `cluster_memberships`,
+            // compute the new centroids, the means of our clusters
             let new_centroids =
                 KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
 
-            // Check convergence condition (very naive, we need an epsilon tolerance here)
+            // Check the convergence condition: if the new centroids,
+            // after the assignment-update cycle, are closer to the old centroids
+            // than a pre-established tolerance we are finished.
             let distance = centroids.sq_l2_dist(&new_centroids).unwrap();
-            has_converged = distance < tolerance;
+            let has_converged = distance < tolerance;
 
             centroids = new_centroids;
+
+            if has_converged {
+                break;
+            }
         }
 
+        // Set `self.centroids` to the outcome of our optimisation process.
         self.centroids = Some(centroids);
     }
 

From cc7de5d7ea2ae5b6fc5c36881c14622fb85dedeb Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:29:25 +0100
Subject: [PATCH 18/22] Add comments

---
 k_means/src/lib.rs | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 33c7f9f..9e6b424 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -88,31 +88,46 @@ impl KMeans {
         self.centroids = Some(centroids);
     }
 
+    /// Given our observations, `X`, and the index to which each observation belongs,
+    /// stored in `cluster_memberships`,
+    /// we want to compute the mean of all observations in each cluster.
     fn compute_centroids<A, B>(
         X: &ArrayBase<A, Ix2>,
         cluster_memberships: &ArrayBase<B, Ix1>,
-        n_clusters: u16,
     ) -> Array2<f64>
     where
         A: Data<Elem = f64>,
         B: Data<Elem = usize>,
     {
         let (_, n_features) = X.dim();
+        // `centroids` is a cluster index -> rolling mean mapping.
+        // We will update it while we iterate over our observations.
         let mut centroids: HashMap<usize, RollingMean> = HashMap::new();
 
+        // We iterate over our observations and cluster memberships in lock-step.
         let iterator = X.genrows().into_iter().zip(cluster_memberships.iter());
         for (row, cluster_index) in iterator {
+            // If we have already encountered an observation that belongs to the
+            // `cluster_index`th cluster, we retrieve the current rolling mean (our new centroid)
+            // and we update it using the current observation.
             if let Some(rolling_mean) = centroids.get_mut(cluster_index) {
                 rolling_mean.accumulate(&row);
             } else {
+                // If we have not yet encountered an observation that belongs to the
+                // `cluster_index`th cluster, we set its centroid to `row`,
+                // initialising our `RollingMean` accumulator.
                 let new_centroid = RollingMean::new(row.to_owned());
+                // .to_owned takes our `row` view as input and returns an owned array.
                 centroids.insert(*cluster_index, new_centroid);
             }
         }
 
-        let mut new_centroids: Array2<f64> = Array2::zeros((n_clusters as usize, n_features));
+        // Convert our `HashMap` into a 2d array.
+        let mut new_centroids: Array2<f64> = Array2::zeros((centroids.len(), n_features));
         for (cluster_index, centroid) in centroids.into_iter() {
             let mut new_centroid = new_centroids.index_axis_mut(Axis(0), cluster_index);
+            // .assign sets each element of `new_centroid`
+            // to the corresponding element in `centroid.current_mean`.
             new_centroid.assign(&centroid.current_mean);
         }
 

From 0f9c482713ee25327694660cd758949cae818305 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:35:53 +0100
Subject: [PATCH 19/22] Refactor to get a simpler function, update rand

---
 k_means/Cargo.toml | 4 ++--
 k_means/src/lib.rs | 8 ++------
 2 files changed, 4 insertions(+), 8 deletions(-)

diff --git a/k_means/Cargo.toml b/k_means/Cargo.toml
index 3a61d25..317f336 100644
--- a/k_means/Cargo.toml
+++ b/k_means/Cargo.toml
@@ -8,5 +8,5 @@ workspace = ".."
 [dependencies]
 ndarray = {version = "0.12"}
 ndarray-stats = {git = "https://github.com/LukeMathWalker/ndarray-stats", branch = "generalized-deviation"}
-ndarray-rand = "0.9"
-rand = "0.6"
+ndarray-rand = "0.10"
+rand = "0.7"
diff --git a/k_means/src/lib.rs b/k_means/src/lib.rs
index 9e6b424..535fd75 100644
--- a/k_means/src/lib.rs
+++ b/k_means/src/lib.rs
@@ -1,7 +1,6 @@
 #![allow(non_snake_case)]
 use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2};
 use ndarray_stats::DeviationExt;
-use rand::distributions::{Distribution, Uniform};
 use std::collections::HashMap;
 
 /// K-means clustering aims to partition a set of observations
@@ -69,7 +68,7 @@ impl KMeans {
             // Update step: using the newly computed `cluster_memberships`,
             // compute the new centroids, the means of our clusters
             let new_centroids =
-                KMeans::compute_centroids(&X, &cluster_memberships, self.n_clusters);
+                KMeans::compute_centroids(&X, &cluster_memberships);
 
             // Check the convergence condition: if the new centroids,
             // after the assignment-update cycle, are closer to the old centroids
@@ -139,11 +138,8 @@ impl KMeans {
         A: Data<Elem = f64>,
     {
         let (n_samples, _) = X.dim();
-        let distribution = Uniform::from(0..n_samples);
         let mut rng = rand::thread_rng();
-        let indices: Vec<_> = (0..n_clusters)
-            .map(|_| distribution.sample(&mut rng))
-            .collect();
+        let indices = rand::seq::index::sample(&mut rng, n_samples, n_clusters as usize).into_vec();
         X.select(Axis(0), &indices)
     }
 

From 20e109f99f41b1539954ab8e4a9eadeb216f1fcc Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:36:51 +0100
Subject: [PATCH 20/22] Avoid deprecation warning

---
 k_means/src/main.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/k_means/src/main.rs b/k_means/src/main.rs
index 500e19c..65e3ba6 100644
--- a/k_means/src/main.rs
+++ b/k_means/src/main.rs
@@ -1,7 +1,7 @@
 #![allow(non_snake_case)]
 use ndarray::{stack, Array, Array2, Axis};
 use ndarray_rand::RandomExt;
-use rand::distributions::Normal;
+use ndarray_rand::rand_distr::Normal;
 
 // Import KMeans from other file ("lib.rs") in this example
 use k_means::KMeans;
@@ -13,8 +13,8 @@ use k_means::KMeans;
 /// looking at the output.
 fn get_data(n_samples: usize, n_features: usize) -> Array2<f64> {
     let shape = (n_samples / 2, n_features);
-    let X1: Array2<f64> = Array::random(shape, Normal::new(1000., 0.1));
-    let X2: Array2<f64> = Array::random(shape, Normal::new(-1000., 0.1));
+    let X1: Array2<f64> = Array::random(shape, Normal::new(1000., 0.1).unwrap());
+    let X2: Array2<f64> = Array::random(shape, Normal::new(-1000., 0.1).unwrap());
     stack(Axis(0), &[X1.view(), X2.view()]).unwrap().to_owned()
 }
 

From 0f962558aecc6f0da27b8fe7148d8b08356b402f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:37:23 +0100
Subject: [PATCH 21/22] Use master branch of ndarray-rand

---
 k_means/Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/k_means/Cargo.toml b/k_means/Cargo.toml
index 317f336..328230c 100644
--- a/k_means/Cargo.toml
+++ b/k_means/Cargo.toml
@@ -7,6 +7,6 @@ workspace = ".."
 
 [dependencies]
 ndarray = {version = "0.12"}
-ndarray-stats = {git = "https://github.com/LukeMathWalker/ndarray-stats", branch = "generalized-deviation"}
+ndarray-stats = {git = "https://github.com/LukeMathWalker/ndarray-stats", branch = "master"}
 ndarray-rand = "0.10"
 rand = "0.7"

From 47351743036727e4aec9f9134cd8cf74e00f1458 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 21 Sep 2019 18:38:10 +0100
Subject: [PATCH 22/22] Add to README

---
 README.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/README.md b/README.md
index e984d89..e97c01f 100644
--- a/README.md
+++ b/README.md
@@ -7,6 +7,7 @@ Each example folder contains a description and instructions on how to run it. Do
 Table of contents:
 
 - [Linear regression](https://github.com/rust-ndarray/ndarray-examples/tree/master/linear_regression)
+- [K-Means clustering](https://github.com/rust-ndarray/ndarray-examples/tree/master/k_means)
 
 ## Contributing