Add absolute tolerance to KMeans algorithm

MassimoCimmino · MassimoCimmino · commit 395501e98e65 · 2025-03-18T09:52:46.000-04:00
diff --git a/IBPSA/Utilities/Clustering/KMeans.mo b/IBPSA/Utilities/Clustering/KMeans.mo
@@ -6,6 +6,7 @@ impure function KMeans "k-means clustering algorithm"
   input Integer n_samples "Number of samples";
   input Integer n_features "Number of features";
   input Real relTol=1e-5 "Relative tolerance on cluster positions";
+  input Real absTol=1e-8 "Absolute tolerance on cluster positions";
   input Integer max_iter=500 "Maximum number of k-means iterations";
   input Integer n_init=10 "Number of runs with randomized centroid seeds";
   input Integer n_cluster_size=0 "Length of the cluster_size output vector";
@@ -16,7 +17,8 @@ impure function KMeans "k-means clustering algorithm"
 protected
   Real old_centroids[n_clusters,n_features] "Previous iteration centroids";
   Real new_centroids[n_clusters,n_features] "Next iteration centroids";
-  Real delta_centroids "Maximum relative displacement of cluster centroids between two k-means iterations";
+  Real relDelta_centroids "Maximum relative displacement of cluster centroids between two k-means iterations";
+  Real absDelta_centroids "Maximum absolute displacement of cluster centroids between two k-means iterations";
   Integer new_labels[n_samples] "Next iteration cluster labels";
   Real new_inertia "Inertia of the samples during the current run";
   Real inertia "Minimum inertia of the samples since first run";
@@ -31,7 +33,8 @@ algorithm
  id := Modelica.Math.Random.Utilities.initializeImpureRandom(seed);
 
   // ---- Perform n_init successive runs of the k-means algorithm
- for run in 1:n_init loop
+  inertia := 0.;
+  for run in 1:n_init loop
     // ---- Select initial centroids at random
     // Select 3 non-repeated data points in the data set
     n := Modelica.Math.Random.Utilities.impureRandomInteger(id,1,n_samples);
@@ -53,8 +56,9 @@ algorithm
 
     // ---- k-means iterations
     k_iter := 0;
-    delta_centroids := 2*relTol;
-    while k_iter < max_iter and delta_centroids > relTol loop
+    relDelta_centroids := 2*relTol;
+    absDelta_centroids := 2*absTol;
+    while k_iter < max_iter and (relDelta_centroids > relTol or absDelta_centroids > absTol) loop
       k_iter := k_iter + 1;
 
       // Find centroid closest to each data point
@@ -71,7 +75,8 @@ algorithm
       end for;
 
       // Re-evaluate position of the centroids
-      delta_centroids := 0;
+      relDelta_centroids := 0;
+      absDelta_centroids := 0;
       for j in 1:n_clusters loop
         n := sum(if new_labels[i]==j then 1 else 0 for i in 1:n_samples);
         new_centroids[j,:] := zeros(n_features);
@@ -84,7 +89,8 @@ algorithm
         else
           new_centroids[j,:] := old_centroids[j,:];
         end if;
-        delta_centroids := max(delta_centroids, sum((new_centroids[j,:] - old_centroids[j,:])./old_centroids[j,:]));
+        relDelta_centroids := max(relDelta_centroids, sum(abs(new_centroids[j,:] - old_centroids[j,:]) ./ (old_centroids[j,:] .+ Modelica.Constants.eps)));
+        absDelta_centroids := max(absDelta_centroids, sum(abs(new_centroids[j,:] - old_centroids[j,:])));
       end for;
       old_centroids := new_centroids;
     end while;
@@ -129,6 +135,13 @@ modifying the constant <code>seed</code>.
 </html>", revisions="<html>
 <ul>
 <li>
+March 18, 2025 by Massimo Cimmino<br/>
+Added absolute tolerance. The algorithm stops when any of the relative and
+absolute tolerances is satisfied. This fixes errors that occur when a centroid
+has a value close to zero on any of its axes. See
+<a href=\"https://github.com/ibpsa/modelica-ibpsa/issues/1985\">#1985</a>. 
+</li>
+<li>
 February 1, 2023, by Michael Wetter:<br/>
 Added <code>impure</code> declaration which is needed for compliance with the Modelica Language Specification,
 and is required by Optimica.
