function [W,class,distances]=kmeans(train,k,p0,p1)
% [W,class]=kmeans(train,k,...)
% [W,class,distances]=kmeans(train,k,...)
%	Perform k-means clustering on the training data, with fixed k size
%	Additional optional flags:
%	  'Graphic'	Show a plot of the error as it evolves
%	  'Verbose'	Print intermediate errors to stdout
%	  'Random'	Initialize the codebook randomly (default: initialize
%			it with the furthest apart input vectors).

Graphic=0;
Verbose=0;
RandomInit=0;
i=0;
while exist(sprintf('p%d',i))==1
  parameter=eval(sprintf('p%d',i));
  if strcmp(lower(parameter),'graphic')
    Graphic=1; 
  elseif strcmp(lower(parameter),'verbose')
    Verbose=1;
  elseif strcmp(lower(parameter),'random')
    RandomInit=1;
  end
  i=i+1; end

if Verbose
  fprintf(1,'k-means clustering, residual:'); 
end

if RandomInit
  % Randomly initialize the codebook by selecting k different vectors from
  % the training dataset
  AllIdx=1:(size(train,2));
  Idx=[];
  for i=0:k-1
    NewIdx=round(rand*(length(AllIdx)-1)+1);
    Idx=[Idx AllIdx(NewIdx)];
    AllIdx=AllIdx([1:NewIdx-1 NewIdx+1:length(AllIdx)]);
  end
  W=train(:,Idx); else
  % Select two most distant vectors as a start
  MaxDist=-Inf;
  for i=1:size(train,2)-1
    for j=i+1:size(train,2)
      Dist=sum((train(:,i)-train(:,j)).^2);
      if Dist>MaxDist
	MaxDist=Dist;
	Idx=[i j];
      end
    end
  end
  % Select the vectors farthest away from already selected centers
  for i=3:k
    MaxDist=-Inf;
    for j=notin(1:size(train,2),Idx)
      Dist=0;
      for k=Idx
	Dist=Dist+sum((train(:,k)-train(:,j)).^2);
      end
      if Dist>MaxDist
	MaxDist=Dist;
	NewIdx=j;
      end
    end
    Idx=[Idx NewIdx];
  end
end
W=train(:,Idx);

% First pass
[dist,class]=mindist(train,W);
errordiff=1;
dist=mean(dist);
prev_dist=dist;
distances=dist;

if Graphic
  % Display the convergence plot
  h=figure;
  hold on
  axis([0 35 0 dist]);
  title('Mean distance from codebook');
  plot(distances);
end
  
% Next passes
while (errordiff > 0),
  % Update the codebook
  for i=1:k,
    cl=~(class-i);
    for j=2:size(W,1),
      cl(j,:)=cl(1,:);
    end
    classVects=train .* cl;
    if sum(cl(1,:)) > 0,
      W(:,i)=sum(classVects')' / sum(cl(1,:));
    end
  end
  
  % Classify etc.
  [dist,class]=mindist(train,W);
  dist=mean(dist);
  errordiff=prev_dist-dist;
  prev_dist=dist;
  distances=[distances dist];
  if Graphic
    plot(distances);
    drawnow;
  end
  if Verbose
    fprintf(1,' %4.2f',dist);
  end
end

if Verbose
  fprintf(1,'\n'); 
end

if Graphic
  close(h)
end

