add medoid silhouette, prepare 0.4.0

Author: kno10

CHANGELOG.md

@@ -2,9 +2,11 @@
 
 For changes to the main Rust package, please see <https://github.com/kno10/rust-kmedoids/blob/main/CHANGELOG.md>
 
-## kmedoids 0.4.0 (2022-09-23)
+## kmedoids 0.4.0 (2022-09-24)
 
-- Add PAMSIL, PAMMEDSIL, FastMSC, FasterMSC (Lars Lenssen)
+- add clustering by optimizing the Silhouette: PAMSIL
+- add medoid silhouette
+- add medoid silhouette clustering: PAMMEDSIL, FastMSC, FasterMSC
 
 ## kmedoids 0.3.3 (2022-07-06)
 

README.md

@@ -3,6 +3,13 @@
 This python package implements k-medoids clustering with PAM and variants of clustering by direct optimization of the (Medoid) Silhouette.
 It can be used with arbitrary dissimilarites, as it requires a dissimilarity matrix as input.
 
+This software package has been introduced in JOSS:
+
+> Erich Schubert and Lars Lenssen  
+> **Fast k-medoids Clustering in Rust and Python**  
+> Journal of Open Source Software 7(75), 4183  
+> <https://doi.org/10.21105/joss.04183> (open access)
+
 For further details on the implemented algorithm FasterPAM, see:
 
 > Erich Schubert, Peter J. Rousseeuw  
@@ -120,10 +127,11 @@ print("Loss with PAM:", pam.loss)
 * PAM (Kaufman and Rousseeuw, 1987) with BUILD and SWAP
 * Alternating optimization (k-means-style algorithm)
 * Silhouette index for evaluation (Rousseeuw, 1987)
-* FasterMSC (Lenssen and Schubert, 2022)
+* **FasterMSC** (Lenssen and Schubert, 2022)
 * FastMSC (Lenssen and Schubert, 2022)
 * PAMSIL (Van der Laan and Pollard, 2003)
 * PAMMEDSIL (Van der Laan and Pollard, 2003)
+* Medoid Silhouette index for evaluation (Van der Laan and Pollard, 2003)
 
 Note that the k-means-like algorithm for k-medoids tends to find much worse solutions.
 

docs/index.rst

@@ -7,6 +7,8 @@ This package is a wrapper around the fast
 `Rust k-medoids package <https://github.com/kno10/rust-kmedoids>`_,
 implementing the FasterPAM and FastPAM algorithms
 along with the baseline k-means-style and PAM algorithms.
+Furthermore, the (Medoid) Silhouette can be optimized by the
+FasterMSC, FastMSC, PAMMEDSIL and PAMSIL algorithms.
 
 All algorithms expect a *distance matrix* and the number of clusters as input.
 They hence can be used with arbitrary dissimilarity functions.
@@ -111,8 +113,10 @@ Implemented Algorithms
 * :ref:`FastMSC<fastmsc>` (Lenssen and Schubert, 2022)
 * :ref:`PAMSIL<pamsil>` (Van der Laan and Pollard, 2003)
 * :ref:`PAMMEDSIL<pammedsil>` (Van der Laan and Pollard, 2003)
+* :ref:`Medoid Silhouette<medoid_silhouette>` (Van der Laan and Pollard, 2003)
 
-Note that the k-means style "alternating" algorithm yields rather poor result quality.
+Note that the k-means style "alternating" algorithm yields rather poor result quality
+(see Schubert and Rousseeuw 2021 for an example and explanation).
 
 .. _FasterPAM:
 
@@ -149,13 +153,6 @@ PAM BUILD
 
 .. autofunction:: pam_build
 
-.. _Silhouette:
-
-Silhouette
-=========
-
-.. autofunction:: silhouette
-
 .. _FasterMSC:
 
 FasterMSC
@@ -184,6 +181,20 @@ PAMMEDSIL
 
 .. autofunction:: pammedsil
 
+.. _Silhouette:
+
+Silhouette
+=========
+
+.. autofunction:: silhouette
+
+.. _MedoidSilhouette:
+
+Medoid Silhouette
+=================
+
+.. autofunction:: medoid_silhouette
+
 .. _KMedoidsResult:
 
 k-Medoids result object
@@ -205,6 +216,13 @@ sklearn-compatible API
 References
 ==========
 
+This software has been published in the Journal of Open-Source Software:
+
+     | Erich Schubert and Lars Lenssen:  
+     | **Fast k-medoids Clustering in Rust and Python**  
+     | Journal of Open Source Software 7(75), 4183  
+     | https://doi.org/10.21105/joss.04183 (open access)
+
 For further details on the implemented algorithm FasterPAM, see:
 
      | Erich Schubert, Peter J. Rousseeuw  
@@ -251,9 +269,3 @@ License: GPL-3 or later
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-FAQ: Why GPL and not Apache/MIT/BSD?
-------------------------------------
-
-Because copyleft software like Linux is what built the open-source community.
-
-Tit for tat: you get to use my code, I get to use your code.

kmedoids/__init__.py

@@ -21,6 +21,11 @@
 
 References:
 
+| Erich Schubert and Lars Lenssen:
+| Fast k-medoids Clustering in Rust and Python
+| Journal of Open Source Software 7(75), 4183
+| https://doi.org/10.21105/joss.04183 (open access)
+
 | Erich Schubert, Peter J. Rousseeuw
 | Fast and Eager k-Medoids Clustering:
 | O(k) Runtime Improvement of the PAM, CLARA, and CLARANS Algorithms
@@ -603,7 +608,7 @@ def silhouette(diss, labels, samples=False, n_cpu=-1):
 	:param n_cpu: number of threads to use (-1: automatic)
 	:type n_cpu: int
 
-	:return: tuple containing the overall silhouette and the individual samples
+	:return: tuple containing the average silhouette and the individual samples
 	:rtype: (float, ndarray)
 	"""
 	import numpy as np, os
@@ -613,7 +618,7 @@ def silhouette(diss, labels, samples=False, n_cpu=-1):
 	if not isinstance(diss, np.ndarray):
 		diss = np.array(diss)
 	if not isinstance(labels, np.ndarray):
-		labels = np.array(labels, dtype=np.uint)
+		labels = np.array(labels, dtype=np.uint64)
 
 	if isinstance(diss, np.ndarray):
 		dtype = diss.dtype
@@ -623,11 +628,11 @@ def silhouette(diss, labels, samples=False, n_cpu=-1):
 		if n_cpu > 1:
 			assert not samples, "samples=true currently may only be used with n_cpu=1"
 			if dtype == np.float32:
-				return _par_silhouette_f32(diss, labels.astype(np.uint64), n_cpu)
+				return (_par_silhouette_f32(diss, labels.astype(np.uint64), n_cpu), [])
 			elif dtype == np.float64:
-				return _par_silhouette_f64(diss, labels.astype(np.uint64), n_cpu)
+				return (_par_silhouette_f64(diss, labels.astype(np.uint64), n_cpu), [])
 			elif dtype == np.int32:
-				return _par_silhouette_i32(diss, labels.astype(np.uint64), n_cpu)
+				return (_par_silhouette_i32(diss, labels.astype(np.uint64), n_cpu), [])
 			elif dtype == np.int64:
 				raise ValueError("Input of int64 is currently not supported, as it could overflow the float64 used internally when computing Silhouette. Use diss.astype(numpy.float64) if that is acceptable and you have the necessary memory for this copy.")
 		else:
@@ -641,6 +646,54 @@ def silhouette(diss, labels, samples=False, n_cpu=-1):
 				raise ValueError("Input of int64 is currently not supported, as it could overflow the float64 used internally when computing Silhouette. Use diss.astype(numpy.float64) if that is acceptable and you have the necessary memory for this copy.")
 	raise ValueError("Input data not supported. Use a numpy array of floats.")
 
+def medoid_silhouette(diss, meds, samples=False):
+	"""Medoid Silhouette index for cluster evaluation.
+
+	The Medoid Silhouette is an approximation to the Silhouette index, that
+	uses the distance to the cluster medoids instead of the average distance
+	to the other cluster members. If every point is assigned to the nearest
+	medoid, the Medoid Silhouette of a point reduces to 1-a/b where a is the
+	distance to the nearest, and b the distance to the second nearest medoid.
+	If b is 0, the Medoid Silhouette is 1.
+
+	This function assumes you already have a distance matrix. It is not necessary
+	to compute a distance matrix to evaluate the medoid silhouette -- only the
+	distances between points and medoids are necessary. If you do not have a
+	distance matrix, simply compute the medoid Silhouette directly, by computing
+	(1) the N x k distance matrix to the medoids, (2) finding the two smallest values
+	for each data point, and (3) computing the average of 1-a/b on these (with 0/0 as 0).
+	This can be implemented in a few lines with numpy easily.
+
+	:param diss: square numpy array of dissimilarities
+	:type diss: ndarray
+	:param meds: medoid indexes (k distinct values in 0 to n-1)
+	:type meds: ndarray of int
+	:param samples: whether to return individual samples or not
+	:type samples: boolean
+
+	:return: tuple containing the average Medoid Silhouette and the individual samples
+	:rtype: (float, ndarray)
+	"""
+	import numpy as np, os
+	from .kmedoids import _medoid_silhouette_i32, _medoid_silhouette_f32, _medoid_silhouette_f64
+
+	if not isinstance(diss, np.ndarray):
+		diss = np.array(diss)
+	if not isinstance(meds, np.ndarray):
+		meds = np.array(meds, dtype=np.uint64)
+
+	if isinstance(diss, np.ndarray):
+		dtype = diss.dtype
+		if dtype == np.float32:
+			return _medoid_silhouette_f32(diss, meds.astype(np.uint64), samples)
+		elif dtype == np.float64:
+			return _medoid_silhouette_f64(diss, meds.astype(np.uint64), samples)
+		elif dtype == np.int32:
+			return _medoid_silhouette_i32(diss, meds.astype(np.uint64), samples)
+		elif dtype == np.int64:
+			raise ValueError("Input of int64 is currently not supported, as it could overflow the float64 used internally when computing Silhouette. Use diss.astype(numpy.float64) if that is acceptable and you have the necessary memory for this copy.")
+	raise ValueError("Input data not supported. Use a numpy array of floats.")
+
 # This is a hack to make sklearn an optional dependency only:
 try:
 	from sklearn.base import BaseEstimator, ClusterMixin, TransformerMixin
@@ -655,8 +708,13 @@ class KMedoids(SKLearnClusterer):
 
 	References:
 
-	| Erich Schubert, Peter J. Rousseeuw
-	| Fast and Eager k-Medoids Clustering:
+	| Erich Schubert and Lars Lenssen:
+	| Fast k-medoids Clustering in Rust and Python
+	| Journal of Open Source Software 7(75), 4183
+	| https://doi.org/10.21105/joss.04183 (open access)
+
+	| Erich Schubert, Peter J. Rousseeuw:
+	| Fast and Eager k-Medoids Clustering
 	| O(k) Runtime Improvement of the PAM, CLARA, and CLARANS Algorithms
 	| Information Systems (101), 2021, 101804
 	| https://doi.org/10.1016/j.is.2021.101804 (open access)

src/lib.rs

@@ -225,6 +225,33 @@ par_silhouette_call!(par_silhouette_f64, f64);
 par_silhouette_call!(par_silhouette_i32, i32);
 // i64 not supported, as the f64 used internally may overflow
 
+macro_rules! medoid_silhouette_call {
+($name:ident, $type: ty) => {
+/// Run the Medoid Silhouette index evaluation for $type precision
+///
+/// :param dist: distance matrix
+/// :type dist: ndarray
+/// :param meds: medoids indexes
+/// :type meds: ndarray
+/// :param samples: return the per-point Medoid Silhouette values
+/// :type samples: bool
+/// :return: Medoid Silhouette evaluation result
+/// :rtype: pair of Medoid Silhouette score and Medoid Silhouette coefficients per point
+#[pyfunction]
+fn $name(dist: PyReadonlyArray2<'_, $type>, meds: PyReadonlyArray1<'_, usize>, samples: bool) -> PyResult<Py<PyAny>> {
+    assert_eq!(dist.ndim(), 2);
+    assert_eq!(dist.shape()[0], dist.shape()[1]);
+    let (sil, sils): (f64, _) = rustkmedoids::medoid_silhouette(&dist.as_array(), &meds.to_vec()?, samples);
+    Python::with_gil(|py| -> PyResult<Py<PyAny>> {
+      Ok((sil, PyArray1::from_vec(py, sils)).to_object(py))
+    })
+}
+}}
+medoid_silhouette_call!(medoid_silhouette_f32, f32);
+medoid_silhouette_call!(medoid_silhouette_f64, f64);
+medoid_silhouette_call!(medoid_silhouette_i32, i32);
+// i64 not supported, as the f64 used internally may overflow
+
 #[pymodule]
 #[allow(unused_variables)]
 fn kmedoids(py: Python, m: &PyModule) -> PyResult<()> {
@@ -263,14 +290,19 @@ fn kmedoids(py: Python, m: &PyModule) -> PyResult<()> {
     m.add("_silhouette_f32", wrap_pyfunction!(silhouette_f32, m)?)?;
     m.add("_silhouette_f64", wrap_pyfunction!(silhouette_f64, m)?)?;
     m.add("_silhouette_i32", wrap_pyfunction!(silhouette_i32, m)?)?;
+    // not supported: m.add("_silhouette_i64", wrap_pyfunction!(silhouette_i64, m)?)?;
     m.add("_par_fasterpam_f32", wrap_pyfunction!(par_fasterpam_f32, m)?)?;
     m.add("_par_fasterpam_f64", wrap_pyfunction!(par_fasterpam_f64, m)?)?;
     m.add("_par_fasterpam_i32", wrap_pyfunction!(par_fasterpam_i32, m)?)?;
     m.add("_par_fasterpam_i64", wrap_pyfunction!(par_fasterpam_i64, m)?)?;
     m.add("_par_silhouette_f32", wrap_pyfunction!(par_silhouette_f32, m)?)?;
     m.add("_par_silhouette_f64", wrap_pyfunction!(par_silhouette_f64, m)?)?;
     m.add("_par_silhouette_i32", wrap_pyfunction!(par_silhouette_i32, m)?)?;
-    // not supported: m.add("_silhouette_i64", wrap_pyfunction!(silhouette_i64, m)?)?;
+    // not supported: m.add("_üarsilhouette_i64", wrap_pyfunction!(par_silhouette_i64, m)?)?;
+    m.add("_medoid_silhouette_f32", wrap_pyfunction!(medoid_silhouette_f32, m)?)?;
+    m.add("_medoid_silhouette_f64", wrap_pyfunction!(medoid_silhouette_f64, m)?)?;
+    m.add("_medoid_silhouette_i32", wrap_pyfunction!(medoid_silhouette_i32, m)?)?;
+    // not supported: m.add("_medoid_silhouette_i64", wrap_pyfunction!(medoid_silhouette_i64, m)?)?;
     Ok(())
 }
 

tests/test_integration.py

@@ -43,14 +43,53 @@ def test_alternating(self):
         assert np.array_equal(alt.medoids, alt_rust[2])
         assert alt.loss == alt_rust[0]
 
+    def test_pamsil(self):
+        dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.float32)
+        pamsil = kmedoids.pamsil(dist, 2)
+        pamsil_rust = kmedoids.kmedoids._pamsil_swap_f32(dist, pamsil.medoids, 100)
+        assert pamsil.loss == 0.5137878787878788
+        assert pamsil.loss == pamsil_rust[0]
+        assert np.array_equal(pamsil.medoids, pamsil_rust[2])
+
+    def test_pammedsil(self):
+        dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.float32)
+        pms = kmedoids.pammedsil(dist, 2)
+        pms_rust = kmedoids.kmedoids._pammedsil_swap_f32(dist, pms.medoids, 100)
+        assert pms.loss == 0.8172727272727272
+        assert pms.loss == pms_rust[0]
+        assert np.array_equal(pms.medoids, pms_rust[2])
+
+    def test_fastmsc(self):
+        dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.float32)
+        fmsc = kmedoids.fastmsc(dist, 2, init="build")
+        fmsc_rust = kmedoids.kmedoids._fastmsc_f32(dist, fmsc.medoids, 100)
+        assert fmsc.loss == 0.8172727272727272
+        assert np.array_equal(fmsc.medoids, fmsc_rust[2])
+        assert fmsc.loss == fmsc_rust[0]
+
+    def test_fastermsc(self):
+        dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.float32)
+        fmsc = kmedoids.fastermsc(dist, 2, init="build")
+        fmsc_rust = kmedoids.kmedoids._fastermsc_f32(dist, fmsc.medoids, 100)
+        assert fmsc.loss == 0.8172727272727272
+        assert np.array_equal(fmsc.medoids, fmsc_rust[2])
+        assert fmsc.loss == fmsc_rust[0]
+
     def test_silhouette(self):
         dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.int32)
         pam = kmedoids.pam(dist, 2)
-        sil = kmedoids.silhouette(dist, pam.labels)
+        sil = kmedoids.silhouette(dist, pam.labels, n_cpu=1)
         par_sil = kmedoids.silhouette(dist, pam.labels, n_cpu=2)
         sil_rust = kmedoids.kmedoids._silhouette_i32(dist, pam.labels, False)
-        assert sil == par_sil
-        assert sil == sil_rust[0]
+        assert sil[0] == par_sil
+        assert sil[0] == sil_rust[0]
+
+    def test_medoid_silhouette(self):
+        dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.float32)
+        pam = kmedoids.fasterpam(dist, 2)
+        sil = kmedoids.medoid_silhouette(dist, pam.labels)
+        sil_rust = kmedoids.kmedoids._medoid_silhouette_f32(dist, pam.labels, False)
+        assert sil[0] == sil_rust[0]
 
     def test_sklearn_interface(self):
         dist = np.array([[0, 2, 3, 4, 5], [2, 0, 6, 7, 8], [3, 6, 0, 9, 10], [4, 7, 9, 0, 11], [5, 8, 10, 11, 0]], dtype=np.int32)