function [X,K,rho,mu,Psi,r,lik] = demo_mog(X,K);

say Mixture of Gaussians Demo
if nargin == 0
  say No dataset given. Would you like to load one or make one (default)?
  cd datasets
  ls *.mat
  cd ..
  dd=input('','s');
  if length(dd)==0
    X = getpoints;
    fname=input('Would you like to save this dataset? If so enter name: ','s');
    if length(fname)>0
      save(['datasets/' fname]);
    end
  else
    load(dd);
  end
  subplot(121);
  plot(X(:,1),X(:,2),'k.','markersize',10); 
  K = input('Number of Gaussian clusters? ')
  dopause = input('Step through demo (yes/no; default yes)? ','s');
  if length(dopause)>0 && dopause(1)=='n'
    numrepeats = 10; % repeat mixture of gaussians EM 10 times.
    dopause = 0;
  else
    numrepeats = 1;
    dopause = 1;
  end
end

m = size(X,1);
n = size(X,2);
if n ~= 2, error('demo only works in 2D'); end

% figure out range of data.
xmin = min(X,[],1);
xmax = max(X,[],1);
xwidth = xmax-xmin;

% colours of clusters
col = colormap;
col = col(ceil((.5:K)/K*size(col,1)),:)*.6;

% keep result of MoG EM for each repeat.
alllik = cell(1,numrepeats); 
allr = cell(1,numrepeats);
allmu = cell(1,numrepeats);
allPsi = cell(1,numrepeats);


for repeat = 1:numrepeats

  % randomly initialize mixture of gaussian
  rho = 1+rand(1,K);
  rho = rho/sum(rho);
  c = randperm(m);
  mu = X(c(1:K),:);
  S = cov(X);
  Psi = S(:,:,ones(1,K))/10;
  r = zeros(m,K);
  lik = zeros(1,40);

  % initial plot
  subplot(121);
  plot(X(:,1),X(:,2),'k.','markersize',10); hold on;
  for c=1:K
    plot(mu(c,1),mu(c,2),'+','markersize',10,'linewidth',3,'color',col(c,:));
  end
  for k=1:K
    plotellipse(mu(k,:),Psi(:,:,k),'linewidth',3,'color',col(k,:));
  end
  hold off
  title(sprintf('MoG EM (repeat %d/%d)',repeat,numrepeats));
  subplot(122);
  cla;
  drawnow
  
  % iterate MoG EM
  for iter = 1:40
    % Run MoG EM for 1 iteration.  
    [rho mu Psi r lik(iter)] = mog(X,K,rho,mu,Psi,r);
  
    % plot 
    subplot(121);
    cla;
    hold on;
    for i = 1:m
      plot(X(i,1),X(i,2),'.','markersize',10,'color',min(r(i,:)*col,1));
    end
    for k=1:K
      plotellipse(mu(k,:),Psi(:,:,k),'linewidth',3,'color',col(k,:));
    end
    hold off;
    title(sprintf('MoG EM (repeat %d/%d)',repeat,numrepeats));
    subplot(122);
    plot(0:iter-1,lik(1:iter),'x-','linewidth',3);
    title('Log likelihood');
    if dopause, pause; else drawnow; end
    % test for convergence
    if iter>=3 && max(max(abs(r-oldr)))<1e-3
      break;
    end
    oldr = r;

  end
  drawnow

  alllik{repeat} = lik(1:iter);
  allr{repeat} = r;
  allmu{repeat} = mu;
  allPsi{repeat} = Psi;
end

% best run of MoG EM
lik = zeros(1,numrepeats);
subplot(122);
cla
hold on
for repeat = 1:numrepeats
  plot(0:length(alllik{repeat})-1,alllik{repeat},'-');
  lik(repeat) = alllik{repeat}(end);
end
hold off
title('MoG EM final log likelihoods');
[bestlik repeat] = max(lik);
r = allr{repeat};
mu = allmu{repeat};
Psi = allPsi{repeat};
subplot(121);
cla;
hold on;
for i = 1:m
  plot(X(i,1),X(i,2),'.','markersize',10,'color',min(r(i,:)*col,1));
end
for k=1:K
  plotellipse(mu(k,:),Psi(:,:,k),'linewidth',3,'color',col(k,:));
end
hold off;
title('Best MoG EM run');

