Repository: phenology/cgc

Issue:
k-means centroids
Centroids returned by the k-means refinement step are averages of the cluster averages. Should one instead return averages computed over all elements of the clusters that are being merged?  

Discussion:

Patch:
diff --git a/cgc/kmeans.py b/cgc/kmeans.py
index 5ba801f..400feb3 100644
--- a/cgc/kmeans.py
+++ b/cgc/kmeans.py
@@ -15,23 +15,23 @@ class KmeansResults(Results):
 
     :var k_value: Optimal K value (value with maximum Silhouette score).
     :type k_value: int
-    :var km_labels: Refined clusters labels. It is a 2D- (for coclustering)
-                    or 3D- (for triclustering) array, with the shape of
-                    `nclusters`. The value at location (band, row, column)
-                    represents the refined cluster label of the corresponding
-                    band/row/column cluster combination.
-    :type km_labels: np.ndarray
+    :var labels: Refined clusters labels. It is a 2D- (for coclustering)
+                 or 3D- (for triclustering) array, with the shape of
+                 `nclusters`. The value at location (band, row, column)
+                 represents the refined cluster label of the corresponding
+                 band/row/column cluster combination.
+    :type labels: np.ndarray
     :var measure_list: List of Silhouette coefficients for all tested k values.
     :type measure_list: np.ndarray
-    :var cl_mean_centroids: Refined cluster averages computed as the centroids
-        of the clusters resulting from the k-means analysis using `k=k_value`.
-        Initally empty clusters are assigned NaN values.
-    :type cl_mean_centroids: np.ndarray
+    :var cluster_averages: Refined cluster averages. They are computed as means
+        over all elements of the co-/tri-clusters assigned to the refined
+        clusters. Initially empty clusters are assigned NaN values.
+    :type cluster_averages: np.ndarray
     """
     k_value = None
-    km_labels = None
+    labels = None
     measure_list = None
-    cl_mean_centroids = None
+    cluster_averages = None
 
 
 class Kmeans(object):
@@ -145,9 +145,8 @@ class Kmeans(object):
         self._compute_statistic_measures()
 
         # Search for value k
-        silhouette_avg_list = np.array(
-            [])  # List of average silhouette measure of each k value
-        kmean_cluster_list = []
+        silhouette_avg_list = np.array([])  # average silhouette measure vs k
+        kmean_label_list = []
         for k in self.k_range:
             # Compute Kmean
             kmean_cluster = KMeans(n_clusters=k,
@@ -157,7 +156,7 @@ class Kmeans(object):
                                               kmean_cluster.labels_)
             silhouette_avg_list = np.append(silhouette_avg_list,
                                             silhouette_avg)
-            kmean_cluster_list.append(kmean_cluster)
+            kmean_label_list.append(kmean_cluster.labels_)
         idx_k = np.argmax(silhouette_avg_list)
         if np.sum(silhouette_avg_list == silhouette_avg_list[idx_k]) > 1:
             idx_k_list = np.argwhere(
@@ -170,31 +169,33 @@ class Kmeans(object):
                     self.k_range[idx_k]))
         self.results.measure_list = silhouette_avg_list
         self.results.k_value = self.k_range[idx_k]
-        self.kmean_cluster = kmean_cluster_list[idx_k]
-        del kmean_cluster_list
+        labels = kmean_label_list[idx_k]
 
-        # Scale back centroids of the "mean" dimension
-        centroids_norm = self.kmean_cluster.cluster_centers_[:, 0]
-        stat_max = np.nanmax(self.stat_measures[:, 0])
-        stat_min = np.nanmin(self.stat_measures[:, 0])
-        mean_centroids = centroids_norm * (stat_max - stat_min) + stat_min
-
-        # Assign centroids to each cluster cell
-        cl_mean_centroids = mean_centroids[self.kmean_cluster.labels_]
-
-        # Reshape the centroids of means to the shape of cluster matrix,
-        # taking into account non-constructive clusters
-        self.results.cl_mean_centroids = np.full(self.nclusters, np.nan)
         indices = np.meshgrid(*[np.unique(cl) for cl in self.clusters],
                               indexing='ij')
-        mask = np.zeros_like(self.results.cl_mean_centroids, dtype=bool)
+        mask = np.zeros(self.nclusters, dtype=bool)
         mask[tuple(indices)] = True
-        self.results.cl_mean_centroids[mask] = cl_mean_centroids
 
         # Make a lookup matrix from un-refined clusters to Kmean clusters
         km_labels = np.full(self.nclusters, np.nan)
-        km_labels[mask] = self.kmean_cluster.labels_
-        self.results.km_labels = km_labels
+        km_labels[mask] = labels
+        self.results.labels = km_labels
+
+        # Calculate the means over the refined clusters
+        cluster_averages = np.full(self.nclusters, np.nan)
+        for label in range(self.results.k_value):
+            label_sum = 0.
+            label_n_elements = 0
+            # Loop over all co-/tri-clusters in the selected refined cluster
+            clusters = np.where(km_labels == label)
+            for cluster in zip(*clusters):
+                idx = [np.where(self.clusters[i] == cluster[i])[0]
+                       for i in range(self.Z.ndim)]
+                label_n_elements += np.prod([len(idx_x) for idx_x in idx])
+                idx = np.meshgrid(*idx, indexing='ij')
+                label_sum += self.Z[tuple(idx)].sum()
+            cluster_averages[clusters] = label_sum / label_n_elements
+        self.results.cluster_averages = cluster_averages
 
         self.results.write(filename=self.output_filename)
         return self.results
@@ -207,12 +208,14 @@ class Kmeans(object):
         """
 
         features = np.zeros((*self.nclusters, 6))
-        features[...,
-                 0] = calculate_cluster_feature(self.Z, np.mean, self.clusters,
-                                                self.nclusters)
-        features[...,
-                 1] = calculate_cluster_feature(self.Z, np.std, self.clusters,
-                                                self.nclusters)
+        features[..., 0] = calculate_cluster_feature(self.Z,
+                                                     np.mean,
+                                                     self.clusters,
+                                                     self.nclusters)
+        features[..., 1] = calculate_cluster_feature(self.Z,
+                                                     np.std,
+                                                     self.clusters,
+                                                     self.nclusters)
         features[..., 2] = calculate_cluster_feature(self.Z,
                                                      np.percentile,
                                                      self.clusters,
@@ -223,19 +226,21 @@ class Kmeans(object):
                                                      self.clusters,
                                                      self.nclusters,
                                                      q=95)
-        features[...,
-                 4] = calculate_cluster_feature(self.Z, np.max, self.clusters,
-                                                self.nclusters)
-        features[...,
-                 5] = calculate_cluster_feature(self.Z, np.min, self.clusters,
-                                                self.nclusters)
+        features[..., 4] = calculate_cluster_feature(self.Z,
+                                                     np.max,
+                                                     self.clusters,
+                                                     self.nclusters)
+        features[..., 5] = calculate_cluster_feature(self.Z,
+                                                     np.min,
+                                                     self.clusters,
+                                                     self.nclusters)
 
-        self.stat_measures = features[~np.isnan(features)].reshape((-1, 6))
+        stat_measures = features[~np.isnan(features)].reshape((-1, 6))
 
         # Normalize all statistics to [0, 1]
-        minimum = self.stat_measures.min(axis=0)
-        maximum = self.stat_measures.max(axis=0)
-        self.stat_measures_norm = np.divide(
-            (self.stat_measures - minimum), (maximum - minimum),
-            out=np.zeros_like((self.stat_measures)),
-            where=(maximum - minimum) != 0)
+        minimum = stat_measures.min(axis=0)
+        maximum = stat_measures.max(axis=0)
+        self.stat_measures_norm = np.divide((stat_measures - minimum),
+                                            (maximum - minimum),
+                                            out=np.zeros_like(stat_measures),
+                                            where=(maximum - minimum) != 0)
diff --git a/docs/kmeans.rst b/docs/kmeans.rst
index 7bfb695..a1983ac 100644
--- a/docs/kmeans.rst
+++ b/docs/kmeans.rst
@@ -65,13 +65,13 @@ The ``compute`` function is then called to run the k-means refinement:
 Results
 -------
 
-The optimal ``k`` value and the k-means-based centroids computed over the co- and tri-cluster averages are stored in the
-KmeansResults object:
+The optimal ``k`` value and the refined cluster averages computed over all elements assigned to the co- and tri-clusters
+are stored in the ``KmeansResults`` object:
 
 .. code-block:: python
 
     results.k_value
-    results.cl_mean_centroids
+    results.cluster_averages
 
 
 .. _k-means clustering: https://en.wikipedia.org/wiki/K-means_clustering


Test patch:
diff --git a/tests/test_kmeans.py b/tests/test_kmeans.py
index 3b6e3dd..492b219 100644
--- a/tests/test_kmeans.py
+++ b/tests/test_kmeans.py
@@ -37,6 +37,29 @@ def init_cocluster_km():
     return km
 
 
+def init_cocluster_km_with_noise():
+    Z = np.array([
+        [0.0, 0.2, 1.0, 1.1],
+        [0.2, 0.0, 1.1, 1.0],
+        [1.0, 1.0, 0.0, 0.0],
+        [1.2, 1.2, 0.1, 0.1],
+        [1.0, 1.2, 0.0, 0.1]
+    ])
+    row_clusters = np.array([0, 0, 1, 1, 1])
+    col_clusters = np.array([0, 0, 1, 1])
+    nrow_clusters, ncol_clusters = 3, 2  # 1 non populated row/col cluster
+    clusters = [row_clusters, col_clusters]
+    nclusters = [nrow_clusters, ncol_clusters]
+    k_range = range(2, 4)
+    kmean_max_iter = 2
+    km = Kmeans(Z=Z,
+                clusters=clusters,
+                nclusters=nclusters,
+                k_range=k_range,
+                kmean_max_iter=kmean_max_iter)
+    return km
+
+
 def init_tricluster_km():
     Z1 = np.array([[0, 0, 1, 1], [0, 0, 1, 1], [1, 1, 0, 0]])
     Z = np.full((4, 3, 4), 0)
@@ -107,14 +130,13 @@ class TestKmeans(unittest.TestCase):
         results = np.array([[0., 0., 0., 0., 0., 0.], [1., 0., 1., 1., 1., 1.],
                             [1., 0., 1., 1., 1., 1.], [0., 0., 0., 0., 0.,
                                                        0.]])
-        self.assertTrue(
-            np.all(results == km.stat_measures))  # First colummn is mean
+        self.assertTrue(np.all(results == km.stat_measures_norm))
 
     def test_kmean_labels_coclustering(self):
         km = init_cocluster_km()
         km.compute()
 
-        cl_km_labels = km.results.km_labels  # Refined labels in clusters
+        cl_km_labels = km.results.labels  # Refined labels in clusters
         row_clusters, col_clusters = np.meshgrid(km.clusters[0],
                                                  km.clusters[1],
                                                  indexing='ij')
@@ -129,13 +151,15 @@ class TestKmeans(unittest.TestCase):
     def test_statistic_triclustering(self):
         km = init_tricluster_km()
         km._compute_statistic_measures()
-        results = np.array([[0., 0., 0., 0., 0., 0.], [1., 0., 1., 1., 1., 1.],
-                            [1., 0., 1., 1., 1., 1.], [0., 0., 0., 0., 0., 0.],
-                            [1., 0., 1., 1., 1., 1.], [2., 0., 2., 2., 2., 2.],
-                            [2., 0., 2., 2., 2., 2.], [1., 0., 1., 1., 1.,
-                                                       1.]])
-        self.assertTrue(
-            np.all(results == km.stat_measures))  # First colummn is mean
+        results = np.array([[0.,   0., 0., 0., 0., 0.],
+                            [0.5, 0., 0.5, 0.5, 0.5, 0.5],
+                            [0.5, 0., 0.5, 0.5, 0.5, 0.5],
+                            [0., 0., 0., 0., 0., 0.],
+                            [0.5, 0., 0.5, 0.5, 0.5, 0.5],
+                            [1., 0., 1., 1., 1., 1.],
+                            [1., 0., 1., 1., 1., 1.],
+                            [0.5, 0., 0.5, 0.5, 0.5, 0.5]])
+        self.assertTrue(np.all(results == km.stat_measures_norm))
 
     def test_kvalues_triclustering(self):
         km = init_tricluster_km()
@@ -146,7 +170,7 @@ class TestKmeans(unittest.TestCase):
         km = init_tricluster_km()
         km.compute()
 
-        cl_km_labels = km.results.km_labels  # Refined labels in clusters
+        cl_km_labels = km.results.labels  # Refined labels in clusters
         band_clusters, row_clusters, col_clusters = np.meshgrid(km.clusters[0],
                                                                 km.clusters[1],
                                                                 km.clusters[2],
@@ -159,15 +183,44 @@ class TestKmeans(unittest.TestCase):
             vk = km.Z[np.where(Z_km_labels == k)]
             self.assertTrue(np.all(vk == vk[0]))
 
-    def test_statistic_centroids_shape(self):
+    def test_centroid_values(self):
         km = init_cocluster_km()
         km.compute()
-        self.assertEqual((3, 2), km.results.cl_mean_centroids.shape)
+        self.assertEqual((3, 2), km.results.cluster_averages.shape)
+        target_centroids = np.array([
+            [0., 1.],
+            [1., 0.],
+            [np.nan, np.nan]
+        ])
+        self.assertTrue(
+            np.allclose(
+                target_centroids,
+                km.results.cluster_averages,
+                equal_nan=True
+            )
+        )
+
+    def test_centroid_values_for_km_with_noise(self):
+        km = init_cocluster_km_with_noise()
+        km.compute()
+        self.assertEqual((3, 2), km.results.cluster_averages.shape)
+        target_centroids = np.array([
+            [0.07, 1.08],
+            [1.08, 0.07],
+            [np.nan, np.nan]
+        ])
+        self.assertTrue(
+            np.allclose(
+                target_centroids,
+                km.results.cluster_averages,
+                equal_nan=True
+            )
+        )
 
     def test_centroids_nan(self):
         km = init_cocluster_km()
         km.compute()
-        self.assertTrue(all(np.isnan(km.results.cl_mean_centroids[2, :])))
+        self.assertTrue(all(np.isnan(km.results.cluster_averages[2, :])))
 
     def test_kvalue_does_not_depend_on_krange_order(self):
         # 4 co-clusters, 2 clusters