Update dynmsc

Author: larslenssen

README.md

@@ -146,10 +146,10 @@ from sklearn.metrics.pairwise import euclidean_distances
 X, _ = fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False)
 X = X[:10000]
 diss = euclidean_distances(X)
-bk = kmedoids.bestk(diss, 100)
-print("Optimal number of clusters according to the Medoid Silhouette:", bk.bestk)
-print("Medoid Silhouette over range of k:", bk.losses)
-print("Range of k:", bk.rangek)
+dm = kmedoids.dynmsc(diss, 100)
+print("Optimal number of clusters according to the Medoid Silhouette:", dm.bestk)
+print("Medoid Silhouette over range of k:", dm.losses)
+print("Range of k:", dm.rangek)
 ```
 
 ### Memory Requirements
@@ -170,7 +170,6 @@ For larger data sets, it is recommended to only cluster a representative sample
 * **FasterMSC** (Lenssen and Schubert, 2022)
 * FastMSC (Lenssen and Schubert, 2022)
 * DynMSC (Lenssen and Schubert, 2023)
-* Bestk (Lenssen and Schubert, 2023)
 * 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)

kmedoids/__init__.py

@@ -111,9 +111,24 @@ def __repr__(self):
 		return f"KMedoidsResult(loss={self.loss}, labels={self.labels}, medoids={self.medoids}, n_iter={self.n_iter}, n_swaps={self.n_swap})"
 
 
-class BestkResult:
+class DynkResult:
 	"""
-	Result of choosing the optimal number of clusters according to the Medoid Silhouette.
+	K-medoids clustering result with automatic number of clusters
+
+	:param loss: Loss of this clustering (sum of deviations)
+	:type loss: float
+
+	:param labels: Cluster assignment
+	:type labels: ndarray
+
+	:param medoids: Chosen medoid indexes
+	:type medoids: ndarray
+
+	:param n_iter: Number of iterations
+	:type n_iter: int
+
+	:param n_swap: Number of swaps performed
+	:type n_swap: int
 
 	:param bestk: Best k by Medoid Silhouette
 	:type bestk: int
@@ -124,13 +139,18 @@ class BestkResult:
 	:param rangek: range of k
 	:type rangek: range
 	"""
-	def __init__(self, bestk, losses, rangek):
+	def __init__(self, loss, labels, medoids, bestk, losses, rangek, n_iter=None, n_swap=None):
+		self.loss = loss
+		self.labels = labels
+		self.medoids = medoids
+		self.n_iter = n_iter
+		self.n_swap = n_swap
 		self.bestk = bestk
 		self.losses = losses
 		self.rangek = rangek
 
 	def __repr__(self):
-		return f"BestkResult(bestk={self.bestk}, losses={self.losses}, rangek={self.rangek})"
+		return f"DynkResult(loss={self.loss}, labels={self.labels}, medoids={self.medoids}, bestk={self.bestk}, losses={self.losses}, rangek={self.rangek}, n_iter={self.n_iter}, n_swaps={self.n_swap})"
 
 def _check_medoids(diss, medoids, init, random_state):
 	"""Check the medoids and random_state parameters."""
@@ -609,8 +629,8 @@ def dynmsc(diss, medoids, max_iter=100, init="random", random_state=None):
 	:param random_state: random seed if no medoids are given
 	:type random_state: int, RandomState instance or None
 
-	:return: k-medoids clustering result
-	:rtype: KMedoidsResult
+	:return: k-medoids clustering with automatic number of clusters
+	:rtype: DynkResult
 	"""
 	import numpy as np
 	from .kmedoids import _dynmsc_f32, _dynmsc_f64
@@ -623,53 +643,9 @@ def dynmsc(diss, medoids, max_iter=100, init="random", random_state=None):
 	if isinstance(diss, np.ndarray):
 		dtype = diss.dtype
 		if dtype == np.float32:
-			return KMedoidsResult(*_dynmsc_f32(diss, medoids.astype(np.uint64), max_iter))
-		elif dtype == np.float64:
-			return KMedoidsResult(*_dynmsc_f64(diss, medoids.astype(np.uint64), max_iter))
-	raise ValueError("Input data not supported. Use a numpy array of floats.")
-
-def bestk(diss, medoids=100, max_iter=100, init="random", random_state=None):
-	"""Optimal number of clusters according to the Medoid Silhouette
-
-	This version uses DynMSC to choose the ptimal number of clusters according
-	to the Medoid Silhouette, that performs DynMSC for k = 2 to the number of input medoids
-	and returns k with the highest Average Medoid Silhouette.
-
-	References:
-
-	| Lars Lenssen, Erich Schubert:
-	| Medoid silhouette clustering with automatic cluster number selection
-	| Information Systems (120), 2024, 102290
-	| <https://doi.org/10.1016/j.is.2023.102290>
-
-	:param diss: square numpy array of dissimilarities
-	:type diss: ndarray
-	:param medoids:  maximum number of clusters to find or existing medoids with length of maximum number of clusters to find
-	:type medoids: int or ndarray
-	:param max_iter: maximum number of iterations
-	:type max_iter: int
-	:param init: initialization method
-	:type init: str, "random", "first" or "build"
-	:param random_state: random seed if no medoids are given
-	:type random_state: int, RandomState instance or None
-
-	:return: Result of choosing the optimal number of clusters according to the Medoid Silhouette
-	:rtype: BestkResult
-	"""
-	import numpy as np
-	from .kmedoids import _bestk_f32, _bestk_f64
-
-	if not isinstance(diss, np.ndarray):
-		diss = np.array(diss)
-
-	medoids = _check_medoids(diss, medoids, init, random_state)
-
-	if isinstance(diss, np.ndarray):
-		dtype = diss.dtype
-		if dtype == np.float32:
-			return BestkResult(*_bestk_f32(diss, medoids.astype(np.uint64), max_iter))
+			return DynkResult(*_dynmsc_f32(diss, medoids.astype(np.uint64), max_iter))
 		elif dtype == np.float64:
-			return BestkResult(*_bestk_f64(diss, medoids.astype(np.uint64), max_iter))
+			return DynkResult(*_dynmsc_f64(diss, medoids.astype(np.uint64), max_iter))
 	raise ValueError("Input data not supported. Use a numpy array of floats.")
 
 def alternating(diss, medoids, max_iter=100, init="random", random_state=None):

src/lib.rs

@@ -181,49 +181,23 @@ macro_rules! dynmsc_call {
 /// :param max_iter: maximum number of iterations
 /// :type max_iter: int
 /// :return: k-medoids clustering result
-/// :rtype: KMedoidsResult
+/// :rtype: DynkResult
 #[pyfunction]
 fn $name(dist: PyReadonlyArray2<'_, $type>, meds: PyReadonlyArray1<'_, usize>, max_iter: usize) -> PyResult<Py<PyAny>> {
     assert_eq!(dist.ndim(), 2);
     assert_eq!(dist.shape()[0], dist.shape()[1]);
     let mut meds = meds.to_vec()?;
-    let (loss, assi, n_iter, n_swap, best_meds, _losses): ($ltype, _, _, _, _, _) = rustkmedoids::dynmsc(&dist.as_array(), &mut meds, max_iter);
+    let maxk = meds.len() + 1;
+    let (loss, assi, n_iter, n_swap, best_meds, losses): ($ltype, _, _, _, _, _) = rustkmedoids::dynmsc(&dist.as_array(), &mut meds, max_iter);
+    let bestk = best_meds.len();
     Python::with_gil(|py| -> PyResult<Py<PyAny>> {
-      Ok((loss, PyArray1::from_vec(py, assi), PyArray1::from_vec(py, best_meds), n_iter, n_swap).to_object(py))
+      Ok((loss, PyArray1::from_vec(py, assi), PyArray1::from_vec(py, best_meds), bestk, PyArray1::from_vec(py, losses), (2..maxk).collect::<Vec<usize>>(), n_iter, n_swap).to_object(py))
     })
 }
 }}
 dynmsc_call!(dynmsc_f32, f32, f64);
 dynmsc_call!(dynmsc_f64, f64, f64);
 
-
-macro_rules! bestk_call {
-($name:ident, $type: ty, $ltype: ty) => {
-/// Run $variant k-medoids clustering function for $type precision
-///
-/// :param dist: distance matrix
-/// :type dist: ndarray
-/// :param meds: initial medoids
-/// :type meds: ndarray
-/// :param max_iter: maximum number of iterations
-/// :type max_iter: int
-/// :return: k-medoids clustering result
-/// :rtype: KMedoidsResult
-#[pyfunction]
-fn $name(dist: PyReadonlyArray2<'_, $type>, meds: PyReadonlyArray1<'_, usize>, max_iter: usize) -> PyResult<Py<PyAny>> {
-    assert_eq!(dist.ndim(), 2);
-    assert_eq!(dist.shape()[0], dist.shape()[1]);
-    let mut meds = meds.to_vec()?;
-    let maxk = meds.len() + 1;
-    let (_loss, _assi, _n_iter, _n_swap, best_meds, losses): ($ltype, _, _, _, _, _) = rustkmedoids::dynmsc(&dist.as_array(), &mut meds, max_iter);
-    Python::with_gil(|py| -> PyResult<Py<PyAny>> {
-      Ok((best_meds.len(), PyArray1::from_vec(py, losses), (2..maxk).collect::<Vec<usize>>()).to_object(py))
-    })
-}
-}}
-bestk_call!(bestk_f32, f32, f64);
-bestk_call!(bestk_f64, f64, f64);
-
 macro_rules! silhouette_call {
 ($name:ident, $type: ty) => {
 /// Run the Silhouette index evaluation for $type precision
@@ -327,8 +301,6 @@ fn kmedoids(py: Python, m: &PyModule) -> PyResult<()> {
     m.add("_fastermsc_f64", wrap_pyfunction!(fastermsc_f64, m)?)?;
     m.add("_dynmsc_f32", wrap_pyfunction!(dynmsc_f32, m)?)?;
     m.add("_dynmsc_f64", wrap_pyfunction!(dynmsc_f64, m)?)?;
-    m.add("_bestk_f32", wrap_pyfunction!(bestk_f32, m)?)?;
-    m.add("_bestk_f64", wrap_pyfunction!(bestk_f64, m)?)?;
     m.add("_pam_swap_f32", wrap_pyfunction!(pam_swap_f32, m)?)?;
     m.add("_pam_swap_f64", wrap_pyfunction!(pam_swap_f64, m)?)?;
     m.add("_pam_swap_i32", wrap_pyfunction!(pam_swap_i32, m)?)?;