function [] = kmeans(fun,varargin)

% KMEANS K-means algorithm exploration tool
%
%   Launch it with KMEANS(DATA,NCLUST) where DATA is the matrix
%   of observations (one observation per row) and NCLUST is the
%   desired number of clusters.
%
%   The clusters are initialized with a heuristic that spreads
%   them randomly around mean(DATA) with standard deviation
%   sqrtm(cov(DATA)).
%
%   If you want to set your own initial clusters, use
%   KMEANS(DATA,MEANS) where MEANS is a cell array containing
%   NCLUST initial mean vectors.
%
%   Example: for two clusters
%     means{1} = [1 2]; means{2} = [3 4];
%     kmeans(data,means);
%


global data mmeans vvars ppriors piConst numClust numPts whereMin whereBefore circle hc hp hb hl cmap;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~ischar(fun),

  data = fun;
  [numPts,dim] = size(data);

  startMean = mean(data);
  startVvars = cov(data);
  startDev = sqrtm(startVvars);

  % Init clusters and likelihood
  if iscell(varargin{1}),

    numClust = length(varargin{1});
    mmeans = varargin{1};
    if (length(varargin)<2),
      for i=1:numClust,
	vvars{i} = startVvars;
      end;
    else,
      vvars = varargin{2};
    end;

    if (length(varargin)<3),
      ppriors(1:numClust) = 1 / numClust;
    else,
      ppriors = varargin{3};
    end;

  else,

    numClust = varargin{1};

    for i=1:numClust,
      mmeans{i} = randn(1,dim) * startDev + startMean;
      vvars{i} = startVvars;
      ppriors(i) = 1 / numClust;
    end;

  end;

  % Make colormap for classification
  cmap = hsv(numClust);

  piConst = -0.5 * numPts * log( (2*pi)^dim );
  % Compute distance
  for i=1:numClust,
    x = data - repmat(mmeans{i},numPts,1);
    dist(i,:) = sum( (x .* x)' );
  end;
  % Find minimum distance
  [minDist,whereMin] = min(dist);
  whereBefore = zeros(size(whereMin));


  % Init figure
  hf = figure('name','Kmeans Explorer','numbertitle','off');

  % Attribution plot
  subplot('position',[ 0 0.05 0.45 0.9 ]);
  %plot(data(:,1),data(:,2),'y+');
  set(gca,'xlim',[0 1000],'ylim',[0 3000],'dataaspectratio',[1 1 1]);
  xlabel('F1 (Hz)'); ylabel('F2 (Hz)');
  title('Attribution of points to clusters');
  grid on; zoom on; hold on;

  circle = [cos(linspace(-pi, pi, 100)') sin(linspace(-pi, pi, 100)')];
  for i=1:numClust,
    subData = data(whereMin==i,:);
    hp(i) = line( subData(:,1), subData(:,2), ...
	'linestyle','none','marker','+','color',cmap(i,:));

    hc(i) = line(mmeans{i}(1),mmeans{i}(2), 10, ...
	'marker','o','markersize',12,'color',[1 1 1],'linewidth',2);
  end;


  % Likelihood plot
  subplot('position',[ 0.51 0.5 0.45 0.45 ]);
  hl = plot(0,NaN,'yo-');
  grid on;
  zoom on;
  set(hl,'markersize',5);
  xlabel('Number of iterations');
  title('Cumulated distance J');


  % Init buttons
  btnWdth = 0.15;
  btnHt = 0.1;

  hb(1) = uicontrol('style','push', ...
      'units','normalized', ...
      'position', [0.5 0.3 btnWdth btnHt], ...
      'backgroundcolor', [0.8 0.8 0.8], ...
      'foregroundcolor', [0 0 0], ...
      'string','Iterate once', ...
      'callback', 'kmeans(''iterate'',1);', ...
      'visible', 'on');

  hb(2) = uicontrol('style','push', ...
      'units','normalized', ...
      'position', [0.65 0.3 btnWdth btnHt], ...
      'backgroundcolor', [0.8 0.8 0.8], ...
      'foregroundcolor', [0 0 0], ...
      'string','Iterate 5 times', ...
      'callback', 'kmeans(''iterate'',5);', ...
      'visible', 'on');

  hb(3) = uicontrol('style','push', ...
      'units','normalized', ...
      'position', [0.8 0.3 btnWdth btnHt], ...
      'backgroundcolor', [0.8 0.8 0.8], ...
      'foregroundcolor', [0 0 0], ...
      'string','Iterate 10 times', ...
      'callback', 'kmeans(''iterate'',10);', ...
      'visible', 'on');

  hb(4) = uicontrol( 'style','text', ...
      'units','normalized', 'position',[0.5 0.15 0.45 0.15], ...
      'string','Push a button to start iterating...', ...
      'foregroundcolor', [1 0.3 0.3], 'backgroundcolor',[0 0 0], ...
      'fontsize',16,'fontweight','bold' );

  hb(5) = uicontrol('style','push', ...
      'units','normalized', ...
      'position', [0.5 0.05 btnWdth btnHt], ...
      'backgroundcolor', [0.8 0.8 0.8], ...
      'foregroundcolor', [0 0 0], ...
      'string','Iterate until conv.', ...
      'callback', 'kmeans(''iterate'',10000);', ...
      'visible', 'on');

  hb(6) = uicontrol('style','push', ...
      'units','normalized', ...
      'position', [0.8 0.05 btnWdth btnHt], ...
      'backgroundcolor', [0.8 0.8 0.8], ...
      'foregroundcolor', [0 0 0], ...
      'string','Close', ...
      'callback', 'kmeans(''stop'');', ...
      'visible', 'on');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(fun,'iterate'),

  nIter = varargin{1};

  set(hb,'enable','off');

  for k=1:nIter,

    % Check for convergence
    if all(all(whereBefore == whereMin)),

      % Send message
      xd = get(hl,'xdata'); yd = get(hl,'ydata');
      str = sprintf('Convergence has been reached in %i iterations.\n\nTotal likelihood at the end is %1.2e.', xd(end), yd(end) );
      set(hb(4),'string',str);

      % Update plots for likelihood
      %for i=1:numClust,
%	% Likelihood
%	subData = data(whereMaxLike==i,:);
%	set( hp(i), 'xdata', subData(:,1), 'ydata',subData(:,2) );
%
%	subplot('position',[ 0 0.05 0.45 0.9 ]);
%	ellip = circle * sqrtm(vvars{i}) + repmat(mmeans{i},100,1);
%	line(ellip(:,1),ellip(:,2),10*ones(size(ellip,1),1), ...
%	    'color',[1 1 1],'linewidth',2);
%	line(mmeans{i}(1),mmeans{i}(2), 10, ...
%	    'marker','+','markersize',10,'color',[1 1 1],'linewidth',2);
%	title('Maximum likelihood classes based on k-means clusters');
%      end;

      drawnow;

      % Send results to base workspace
      assignin('base','kmeans_result_means',mmeans);
      assignin('base','kmeans_result_vars',vvars);
      assignin('base','kmeans_result_priors',ppriors);
      disp(' ');
      disp([ 'KMEANS: resulting means, variances and priors are now stored in' ...
	    ' the workspace variables kmeans_result_means,' ...
	    ' kmeans_result_vars and kmeans_result_priors.' ]);

      break;

    else,

      % Update clusters and plots for distance
      for i=1:numClust,
	% Distance
	subData = data(whereMin==i,:);
	set( hp(i), 'xdata', subData(:,1), 'ydata',subData(:,2) );

	% Update clusters
	mmeans{i} = mean( subData );
	vvars{i} = cov( subData );
	ppriors(i) = size(subData,1) / numPts;

	set( hc(i), 'xdata', mmeans{i}(1), 'ydata', mmeans{i}(2) );

      end;

      whereBefore = whereMin;
      for i=1:numClust,
	% Compute distance
	x = data - repmat(mmeans{i},numPts,1);
	dist(i,:) = sum( (x .* x)' );
	%% Compute log-likelihood
	%invSig = inv(vvars{i});
	%logLike(i,:) = -0.5 * ( sum( ((x*invSig).*x)' ) );
	%logLike(i,:) =  logLike(i,:) - 0.5 * log(det(vvars{i})) + log(ppriors(i));
      end;
      % Find minimum distance
      [minDist,whereMin] = min(dist);
      %% Find max likelihood
      %[maxLike,whereMaxLike] = max(logLike);

      % Update likelihood plot
      xd = get(hl,'xdata'); xd = [xd  xd(end)+1];
      %yd = get(hl,'ydata'); yd = [yd (sum(maxLike) + piConst)];
      yd = get(hl,'ydata'); yd = [yd sum(minDist)];
      set(hl,'xdata',xd,'ydata',yd);

      drawnow;

    end;

  end;

  if all(all(whereBefore == whereMin)),
    set(hb([4 6]),'enable','on');
  else,
    set(hb,'enable','on');
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(fun,'stop'),

  clear data mmeans vvars ppriors piConst numClust numPts whereMin whereBefore circle hc hp hb hl cmap;
  close(gcf);

else,

  error('Unknown function.');

end;
