--- kmeans.m	2015-01-09 14:21:24.715376108 +1100
+++ kmeans.m	2015-07-02 10:55:22.483905178 +1000
@@ -1,5 +1,6 @@
 ## Copyright (C) 2011 Soren Hauberg <soren@hauberg.org>
 ## Copyright (C) 2012 Daniel Ward <dwa012@gmail.com>
+## Copyright (C) 2015 Lachlan Andrew <lachlanbis@gmail.com>
 ##
 ## This program is free software; you can redistribute it and/or modify it under
 ## the terms of the GNU General Public License as published by the Free Software
@@ -15,9 +16,11 @@
 ## this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{idx}, @var{centers}] =} kmeans (@var{data}, @var{k}, @var{param1}, @var{value1}, @dots{})
+## @deftypefn {Function File} {[@var{idx}, @var{centers}, @var{sumd}, @var{D}] =} kmeans (@var{data}, @var{k}, @var{param1}, @var{value1}, @dots{})
 ## K-means clustering.
 ##
+## Accepts the same values for parameters 'start', 'replicates', 'MaxIter',
+## 'distance' and 'emptyaction' as Matlab.
 ## @seealso{linkage}
 ## @end deftypefn
 
@@ -26,102 +29,261 @@
 
   ## defaults for options
   emptyaction = "error";
-  start       = "sample";
+  start       = "plus";
+  replicates  = 1;
+  MaxIter     = 100;
+  distance    = "sqeuclidean";
+
+  ReplicatesSetExplicitly = false;
+
+  ## Remove rows containing NaN / NA, but record which rows are used
+  data_idx = !any (isnan (data), 2);
+  original_rows = rows (data);
+  data = data(data_idx,:);
 
   #used for getting the number of samples
   nRows = rows (data);
 
-  ## used to hold the distances from each sample to each class
-  D = zeros (nRows, k);
-
   #used for convergence of the centroids
   err = 1;
 
-  #initial sum of distances
-  sumd = Inf;
-
-  ## Input checking, validate the matrix and k
+  ## Input checking, validate the matrix
   if (!isnumeric (data) || !ismatrix (data) || !isreal (data))
     error ("kmeans: first input argument must be a DxN real data matrix");
-  elseif (!isscalar (k))
-    error ("kmeans: second input argument must be a scalar");
+  elseif (!isnumeric (k))
+    error ("kmeans: second argument must be numeric");
   endif
 
-  if (length (varargin) > 0)
-    ## check for the 'emptyaction' property
-    found = find (strcmpi (prop, "emptyaction") == 1);
-    switch (lower (prop{found+1}))
-      case "singleton"
-        emptyaction = "singleton";
+  ## Parse options
+  while (length (prop) > 0)
+    if (length (prop) < 2)
+      error ("kmeans: Option '%s' has no argument", prop{1});
+    endif
+    switch (lower (prop{1}))
+      case "emptyaction" emptyaction = prop{2};
+      case "start"       start       = prop{2};
+      case "maxiter"     MaxIter     = prop{2};
+      case "distance"    distance    = prop{2};
+
+      case "replicates"  replicates  = prop{2};
+                         ReplicatesSetExplicitly = true;
+
+      case {"display", "onlinephase", "options"}
+        warning ("kemans: Ignoring unimplemented option '%s'", prop{1});
+
       otherwise
-        error ("kmeans: unsupported empty cluster action parameter");
+        error ("kmeans: Unknown option %s", prop{1});
     endswitch
-  endif
+    prop = {prop{3:end}};
+  endwhile
+
+  ## Process options
+
+  ## check for the 'emptyaction' property
+  switch (emptyaction)
+    case {"singleton", "error", "drop"}
+      ;
+    otherwise
+      error ("kmeans: unsupported empty cluster action parameter '%s'", emptyaction);
+  endswitch
+
+  ## check for the 'replicates' property
+  try
+    if isempty (1:replicates)
+      error ("kmeans: invalid number of replicates");
+    endif
+  catch
+    error ("kmeans: invalid number of replicates");
+  end_try_catch
+
+  ## check for the 'MaxIter' property
+  try
+    if (1 < MaxIter)
+    endif
+  catch
+    error ("kmeans: invalid MaxIter");
+  end_try_catch
 
   ## check for the 'start' property
   switch (lower (start))
-    case "sample"
-      idx = randperm (nRows) (1:k);
-      centers = data (idx, :);
+    case {"sample", "plus", "cluster"}
+      start = lower (start);
+    case {"uniform"}
+      start = "uniform"
+      min_data = min (data);
+      range = max (data) - min_data;
     otherwise
-      error ("kmeans: unsupported initial clustering parameter");
+      if (!ismatrix (start) || !isnumeric (start))
+        error ("kmeans: invalid start parameter");
+      endif
+      if (isempty (k))
+        k = rows (start)
+      elseif (rows (start) != k)
+        error ("kmeans: Number of initializers (%d) should match number of centroids (%d)", rows (start), k);
+      endif
+      if (ReplicatesSetExplicitly)
+        if (replicates != size (start, 3))
+           error ("kmeans: The third dimension of the initializer (%d) should match the number of replicates (%d)", size (start, 3), replicates);
+        endif
+      else
+        replicates = size (start, 3);
+      endif
   endswitch
 
-  ## Run the algorithm
-  while err > .001
-    ## Compute distances
-    for i = 1:k
-      D (:, i) = sumsq (data - repmat (centers(i, :), nRows, 1), 2);
-    endfor
+  ## check for the 'distance' property
+  ## dist  returns the distance btwn each row of matrix x and a row vector c
+  switch (lower (distance))
+    case "sqeuclidean"
+      dist = @(x, c) (sumsq (x - repmat (c, rows (x), 1), 2));
+      centroid  = @(x) (mean (x,1));
+    case "cityblock"
+      dist = @(x, c) (sum (abs (x - repmat (c, rows (x), 1)), 2));
+      centroid  = @(x) (median (x,1));
+    case "cosine"
+        ## Pre-normalize all data.
+        ## (when Octave implements normr, will use  data = normr (data) )
+      for i = 1:rows (data)
+        data(i,:) = data(i,:) / sqrt (sumsq (data(i,:)));
+      endfor
+      dist = @(x, c) (1 - (x * c') ./ sqrt (sumsq (c)));
+      centroid = @(x) (mean (x,1));   ## already normalized
+    case "correlation"
+        ## Pre-normalize all data.
+      data = data - mean (data, 2);
+        ## (when Octave implements normr, will use  data = normr (data) )
+      for i = 1:rows (data)
+        data(i,:) = data(i,:) / sqrt (sumsq (data(i,:)));
+      endfor
+
+      dist = @(x, c) (1 - (x * (c-mean (c))') ./ sqrt (sumsq (c-mean (c))));
+      centroid  = @(x) (mean (x,1));   ## already normalized
+    case "hamming"
+      dist = @(x, c) (sum (x != repmat (c, rows (x), 1), 2));
+      centroid  = @(x) (median (x,1));
+    otherwise
+      error ("kmeans: unsupported distance parameter");
+  endswitch
+
+  ## Done processing options
+  ########################################
+
+  ## Now that  k  has been set (possibly by 'replicates' option), check/use it.
+
+  if (!isscalar (k))
+    error ("kmeans: second input argument must be a scalar");
+  endif
 
-    ## Classify
-    [tmp, classes] = min (D, [], 2);
+  ## used to hold the distances from each sample to each class
+  D = zeros (nRows, k);
 
-    ## Calculate new centroids
+  best = Inf;
+  best_centers = [];
+  for rep = 1:replicates
+    ## check for the 'start' property
+    switch (lower (start))
+      case "sample"
+        idx = randperm (nRows, k);
+        centers = data(idx, :);
+      case "plus"                       # k-means++, by Arthur and Vassilios(?)
+        centers(1,:) = data(randi (nRows),:);
+	d = inf (nRows, 1);       # Distance to nearest centroid so far
+	for i = 2:k
+	  d = min (d, dist (data, centers(i-1, :)));
+	  centers(i,:) = data(find (cumsum (d) > rand * sum (d), 1), :);
+	endfor
+      case "cluster"
+        idx = randperm (nRows, ceil (nRows/10));
+        [~, centers] = kmeans (data(idx), k, "start", "sample", "distance", distance);
+      case "uniform"
+	# vectorised 'min_data + range .* rand'
+        centers = bsxfun(@plus, min_data, bsxfun(@times, range, rand (k, columns (data))));
+      otherwise
+        centers = start(:,:,rep);
+    endswitch
+
+    ## Run the algorithm
+    iter = 1;
+
+        ## Classify once before the loop; to set sumd, and  if  MaxIter == 0
+    ## Compute distances and classify
     for i = 1:k
-      ## Get binary vector indicating membership in cluster i
-      membership = (classes == i);
-      ## Check for empty clusters
-      if (sum (membership) == 0)
-        switch emptyaction
-          ## if 'singleton', then find the point that is the
-          ## farthest and add it to the empty cluster
-          case 'singleton'
-           idx=maxCostSampleIndex (data, centers(i,:));
-           classes(idx) = i;
-           membership(idx)=1;
-         ## if 'error' then throw the error
-          otherwise
-            error ("kmeans: empty cluster created");
-        endswitch
-     endif ## end check for empty clusters
-
-      ## update the centroids
-      members = data(membership, :);
-      centers(i, :) = sum(members,1)/size(members,1);
+      D (:, i) = dist (data, centers(i, :));
     endfor
+    [~, classes] = min (D, [], 2);
+    sumd = objCost (D, classes);
 
-    ## calculate the difference in the sum of distances
-    err  = sumd - objCost (data, classes, centers);
-    ## update the current sum of distances
-    sumd = objCost (data, classes, centers);
-  endwhile
-endfunction
+    while (err > .001 && iter++ <= MaxIter)
+      ## Calculate new centroids
+      replacedCentroids = [];        ## Used by "emptyaction = singleton"
+      for i = 1:k
+        ## Get binary vector indicating membership in cluster i
+        membership = (classes == i);
+
+        ## Check for empty clusters
+        if (!any (membership))
+          switch emptyaction
+            ## if 'singleton', then find the point that is the
+            ## farthest from any centroid (and not replacing an empty cluster
+            ## from earlier in this pass) and add it to the empty cluster
+            case 'singleton'
+             available = setdiff(1:nRows, replacedCentroids);
+             [~, idx] = max (min (D(available,:)'));
+             idx = available(idx);
+             replacedCentroids = [replacedCentroids, idx];
+
+             classes(idx) = i;
+             membership(idx)=1;
+
+           ## if 'drop' then set C and D to NA
+           case 'drop'
+            centers(i,:) = NA;
+            D(i,:) = NA;
+
+           ## if 'error' then throw the error
+            otherwise
+              error ("kmeans: empty cluster created");
+          endswitch
+       endif ## end check for empty clusters
+
+        ## update the centroids
+        if (any (membership))      ## if we didn't "drop" the cluster
+          centers(i, :) = centroid (data(membership, :));
+        endif
+      endfor
+
+      ## Compute distances
+      for i = 1:k
+        D (:, i) = dist (data, centers(i, :));
+      endfor
+
+      ## Classify
+      [~, classes] = min (D, [], 2);
+
+      ## calculate the difference in the sum of distances
+      new_sumd = objCost (D, classes);
+      err  = sum (sumd - new_sumd);
+      ## update the current sum of distances
+      sumd = new_sumd;
+    endwhile
+    if (sum (sumd) < sum (best) || isinf (best))
+      best = sumd;
+      best_centers = centers;
+    endif
+  endfor
+  centers = best_centers;
+  sumd = best;
 
-## calculate the sum of distances
-function obj = objCost (data, classes, centers)
-  obj = 0;
-    for i=1:rows (data)
-      obj = obj + sumsq (data(i,:) - centers(classes(i),:));
-    endfor
+  final_classes = NA (original_rows,1);
+  final_classes(data_idx) = classes;        ## other positions already NaN / NA
+  classes = final_classes;
 endfunction
 
-function idx = maxCostSampleIndex (data, centers)
-  cost = 0;
-  for idx = 1:rows (data)
-    if cost < sumsq (data(idx,:) - centers)
-      cost = sumsq (data(idx,:) - centers);
-    endif
+## calculate the sum of within-class distances
+function obj = objCost (D, classes)
+  obj = zeros (1,columns (D));
+  for i = 1:columns (D)
+    idx = (classes == i);
+    obj(i) = sum (D(idx,i));
   endfor
 endfunction
 
