Add a minimum k to DynMSC

Author: larslenssen

README.md

@@ -141,7 +141,7 @@ print("Loss with PAM:", pam.loss)
 ### Choose the optimal number of clusters
 
 This package includes DynMSC, an algorithm that optimizes the Medoid Silhouette,
-and chooses the "optimal" number of clusters in a range of 2..kmax.
+and chooses the "optimal" number of clusters in a range of kmin..kmax.
 Beware that if you allow a too large kmax, the optimum result will likely have many
 one-elemental clusters. A too high kmax may mask more desirable results, hence it
 is recommended that you choose only 2-3 times the number of clusters you expect as maximum.
@@ -153,7 +153,9 @@ 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)
-dm = kmedoids.dynmsc(diss, 100)
+kmin = 10
+kmax = 20
+dm = kmedoids.dynmsc(diss, kmax, kmin)
 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)

docs/index.rst

@@ -118,7 +118,7 @@ Choosing the optimal number of clusters
 ---------------------------------------
 
 This package includes :ref:`DynMSC<dynmsc>`, an algorithm that optimizes the Medoid Silhouette,
-and chooses the "optimal" number of clusters in a range of 2..kmax.
+and chooses the "optimal" number of clusters in a range of kmin..kmax.
 Beware that if you allow a too large kmax, the optimum result will likely have many
 one-elemental clusters. A too high kmax may mask more desirable results, hence it
 is recommended that you choose only 2-3 times the number of clusters you expect as maximum.
@@ -131,7 +131,9 @@ is recommended that you choose only 2-3 times the number of clusters you expect
     X, _ = fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False)
     X = X[:10000]
     diss = euclidean_distances(X)
-    dm = kmedoids.dynmsc(diss, 100)
+    kmin = 10
+    kmax = 20
+    dm = kmedoids.dynmsc(diss, kmax, kmin)
     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)

kmedoids/__init__.py

@@ -613,11 +613,11 @@ def fastermsc(diss, medoids, max_iter=100, init="random", random_state=None):
 			return KMedoidsResult(*_fastermsc_f64(diss, medoids.astype(np.uint64), max_iter))
 	raise ValueError("Input data not supported. Use a numpy array of floats.")
 
-def dynmsc(diss, medoids, max_iter=100, init="random", random_state=None):
+def dynmsc(diss, medoids, minimum_k=2, max_iter=100, init="random", random_state=None):
 	"""DynMSC clustering
 
 	This is a version of FasterMSC with automatic cluster number selection, that
-	performs FasterMSC for k = 2 to the number of input medoids and returns
+	performs FasterMSC for a minimum k to the number of input medoids and returns
 	the clustering with the highest Average Medoid Silhouette.
 
 	References:
@@ -636,6 +636,8 @@ def dynmsc(diss, medoids, max_iter=100, init="random", random_state=None):
 	:type max_iter: int
 	:param init: initialization method
 	:type init: str, "random", "first" or "build"
+	:param minimum_k: minimum number of clusters to find
+	:type minimum_k: int
 	:param random_state: random seed if no medoids are given
 	:type random_state: int, RandomState instance or None
 
@@ -650,12 +652,14 @@ def dynmsc(diss, medoids, max_iter=100, init="random", random_state=None):
 
 	medoids = _check_medoids(diss, medoids, init, random_state)
 
+	if medoids.shape[0] < minimum_k:
+		raise ValueError("Maximum k should be at least minimum k.")
 	if isinstance(diss, np.ndarray):
 		dtype = diss.dtype
 		if dtype == np.float32:
-			return DynkResult(*_dynmsc_f32(diss, medoids.astype(np.uint64), max_iter))
+			return DynkResult(*_dynmsc_f32(diss, medoids.astype(np.uint64), minimum_k, max_iter))
 		elif dtype == np.float64:
-			return DynkResult(*_dynmsc_f64(diss, medoids.astype(np.uint64), max_iter))
+			return DynkResult(*_dynmsc_f64(diss, medoids.astype(np.uint64), minimum_k, 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):
@@ -933,7 +937,7 @@ def fit(self, X, y=None):
 		elif self.method == "fastermsc":
 			result = fastermsc(X, self.n_clusters, self.max_iter, self.init, random_state=self.random_state)
 		elif self.method == "dynmsc":
-			result = dynmsc(X, self.n_clusters, self.max_iter, self.init, random_state=self.random_state)
+			result = dynmsc(X, self.n_clusters, 2, self.max_iter, self.init, random_state=self.random_state)
 		elif self.method == "fastmsc":
 			result = fastmsc(X, self.n_clusters, self.max_iter, self.init, random_state=self.random_state)
 		elif self.method == "pamsil":

src/lib.rs

@@ -183,15 +183,15 @@ macro_rules! dynmsc_call {
 /// :return: k-medoids clustering result
 /// :rtype: DynkResult
 #[pyfunction]
-fn $name(dist: PyReadonlyArray2<'_, $type>, meds: PyReadonlyArray1<'_, usize>, max_iter: usize) -> PyResult<Py<PyAny>> {
+fn $name(dist: PyReadonlyArray2<'_, $type>, meds: PyReadonlyArray1<'_, usize>, minimum_k: 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);
+    let (loss, assi, n_iter, n_swap, best_meds, losses): ($ltype, _, _, _, _, _) = rustkmedoids::dynmsc(&dist.as_array(), &mut meds, minimum_k, 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), bestk, PyArray1::from_vec(py, losses), (2..maxk).collect::<Vec<usize>>(), 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), (minimum_k..maxk).collect::<Vec<usize>>(), n_iter, n_swap).to_object(py))
     })
 }
 }}