function [ c, ptoc, v_data ] = km ( dim, n, p, k )

%*****************************************************************************80
%
%% KM carries out the K-Means clustering algorithm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM, the dimension of the data.
%
%    Input, integer N, the number of data objects.
%
%    Input, real P(DIM,N), the data.
%
%    Input, integer K, the number of clusters to be formed.
%
%    Output, real C(DIM,K), the cluster centers.
%
%    Output, integer PTOC(N), indicates the cluster (between 1 and K)
%    to which each data object was assigned.
%
%    Output, real V_DATA(*), a list of the variances that were computed
%    at each stage of the iteration.  This data can be plotted by the
%    command "plot ( v_data )".
%

%
%  Initialize the cluster centers.
%
  c = kmeans_initialize ( dim, n, p, k );
%
%  Repeatedly update clusters and centers til no change.
%
  v = -1;
  v_data = [];

  while ( 1 )

    ptoc = kmeans_update_clusters ( dim, n, p, k, c );

    c = kmeans_update_centers ( dim, n, p, k, ptoc );

    v_old = v;
    v = kmeans_variance ( dim, n, p, k, c, ptoc );

    v_data = [ v_data, v ];

    if ( v == v_old )
      break
    end

  end

  return
end
function c = kmeans_initialize ( dim, n, p, k )

%*****************************************************************************80
%
%% KMEANS_INITIALIZE randomly chooses K data values for cluster centers.
%
%  Get a random permutation of the integers 1 through N.
%
  perm = randperm ( n );
%
%  Use the first K entries of the permutation to select K 
%  data values to use as cluster centers.
%
  c(1:dim,1:k) = p(1:dim,perm(1:k));

  return
end
function c = kmeans_update_centers ( dim, n, p, k, ptoc )

%*****************************************************************************80
%
%% KMEANS_UPDATE_CENTERS resets the cluster centers to the data averages.
%
  for j = 1 : k
    index = find ( ptoc(1:n) == j );
    nj = length ( index );
    c(1:dim,j) = sum ( p(1:dim,index), 2 ) / nj;
  end

  return
end
function ptoc = kmeans_update_clusters ( dim, n, p, k, c )

%*****************************************************************************80
%
%% KMEANS_UPDATE_CLUSTERS assigns data to clusters based on the centers.
%
  for i = 1 : n
    for j = 1 : k
      dist(j) = norm ( p(1:dim,i) - c(1:dim,j) );
    end
    [ dummy, ptoc(i) ] = min ( dist(1:k) );
  end

  return
end
function v = kmeans_variance ( dim, n, p, k, c, ptoc )

%*****************************************************************************80
%
%% KMEANS_VARIANCE computes the variance of the K-means clustering.
%
  v = 0.0;
  for i = 1 : n
    j = ptoc(i);
    v = v + ( norm ( p(1:dim,i) - c(1:dim,j) ) )^2;
  end

  return
end
