# HG changeset patch
# User cuantiar <info@cuantiar.com>
# Date 1462631952 18000
#      Sat May 07 09:39:12 2016 -0500
# Node ID 8d3af0f9cc51093816f828f212d22cbf01a46cb8
# Parent  42ab417722d3b59f6844284ac98cb073dc916075
make Octave the coolest software ever

diff -r 42ab417722d3 -r 8d3af0f9cc51 inst/kmeans.m
--- a/inst/kmeans.m	Thu Apr 14 15:46:32 2016 +0100
+++ b/inst/kmeans.m	Sat May 07 09:39:12 2016 -0500
@@ -1,3 +1,12 @@
+
+
+    ### modifications by "cuantiar" 26 March 2016
+    ### 3 modifications in total
+    ### modifications No.1, No.2 in lines: 117-139
+    ### modifications No.3 in lines: 145-181
+
+
+
 ## Copyright (C) 2011 Soren Hauberg <soren@hauberg.org>
 ## Copyright (C) 2012 Daniel Ward <dwa012@gmail.com>
 ##
@@ -63,7 +72,7 @@
   ## check for the 'start' property
   switch (lower (start))
     case "sample"
-      idx = randperm (nRows) (1:k);
+      idx = randperm (nRows) (1:k); %% notes
       centers = data (idx, :);
     otherwise
       error ("kmeans: unsupported initial clustering parameter");
@@ -104,19 +113,72 @@
     endfor
 
     ## calculate the difference in the sum of distances
-    err  = sumd - objCost (data, classes, centers);
+
+    ###########################################################################
+    ### modification by "cuantiar" 26 March 2016
+    ### 2 modifications in here
+    ###########################################################################
+
+    ### modification No.1:
+    ### "err" represents a derivative and it may be better to terminate the iterations
+    ### when "err ~ 0" but not when "err <~ 0" (as it was in the code). 
+    ### as it was "err" can take negative values of large absolute value
+    ### and terminate the convergence prematurely
+    ### in the modifications have replaced "err" by "abs(err)"
+
+    ### modification No.2:
+    ### it does not need to calculate the function "objCost" twice, 
+    ### but only once and store its value in "objCostValue" to be reused in the following line
+
+    objCostValue    =  objCost (data, classes, centers);   
+    %%%% err  = sumd - objCost (data, classes, centers);
+    err             = abs(sumd - objCostValue);
     ## update the current sum of distances
-    sumd = objCost (data, classes, centers);
+    ### sumd = objCost (data, classes, centers);
+    sumd            = objCostValue;
+    ###########################################################################
   endwhile
 endfunction
 
 ## calculate the sum of distances
+
+    ###########################################################################
+    ### modification by "cuantiar" 26 March 2016
+    ### 1 modifications in here
+    ###########################################################################
+
+    ### modification No.3:
+
+###########################################################################
+### modification by "cuantiar" 26 March 2016
+### function replaced by the one below
+###########################################################################
+%function obj = objCost (data, classes, centers)
+%  obj = 0;
+%    for i=1:rows (data)
+%      obj = obj + sumsq (data(i,:) - centers(classes(i),:));
+%    endfor
+%endfunction
+
+    ### the "objCost" written in a more vector form, to be faster 
+    ### the non-vector loop is written for the number of classes 
+    ### instead for the number of data-points (as in the previous function version)
+    ### (since the number of classes is much less than the number of points)
+    ### in the test-case attached, the modified code was about x 121 faster (that depends on the data-set)
+#############################
+### objCost
+## calculate the sum of distances
+## over all classes and all data-points
+#############################
 function obj = objCost (data, classes, centers)
-  obj = 0;
-    for i=1:rows (data)
-      obj = obj + sumsq (data(i,:) - centers(classes(i),:));
-    endfor
+    obj = 0.0;
+    for k = 1: size(centers, 1)
+        ck  = find(classes == k);
+        ds  = sumsq(data(ck, :) - centers(k, :))
+        obj = obj + sum(ds);
+    end %for
 endfunction
+#############################
 
 function idx = maxCostSampleIndex (data, centers)
   cost = 0;
@@ -129,8 +191,8 @@
 
 %!demo
 %! ## Generate a two-cluster problem
-%! C1 = randn (100, 2) + 1;
-%! C2 = randn (100, 2) - 1;
+%! C1 = randn (10000, 2) + 1;
+%! C2 = randn (10000, 2) - 1;
 %! data = [C1; C2];
 %!
 %! ## Perform clustering
