function lyngby_km_plot_scat(Centers, varargin)

% lyngby_km_plot_scat  - K-means scattergram of cluster centers
%
%       function lyngby_km_plot_scat(Centers, ...
%          'PropertyName', 'PropertyValue')
% 
%       Input:    Centers  Matrix with cluster centers
%
%       Property: Assign   Vector with assignment corresponding to the
%                          row index of the 'Data' datamatrix
%                 Data     Original data, the number of columns should
%                          equal that in Centers 
%                 Labels   Label (eg. known labels) that should
%                          correspond to the rows 'Data' datamatrix.
%
%       This function plots a scattergram of the cluster centers from
%       a K-means clustering analysis. The cluster center matrix is
%       singular value decomposed and the centers are projected to a
%       2D subspace spanned by the principal components. 
%
%       If 'Data' are given (the original data matrix) then these are
%       projected to the subspace. If 'Assign' is given then the
%       projected elements in 'Data' are colored according to their
%       cluster membership. Furthermore if 'Labels' are given the
%       elements in 'Data' are rendered with different points. 
%
%       Example: 
%         load iris.txt
%         [centers, assign] = lyngby_km_main(iris, 'type', 'mean', 
%            'standardization', 'std', 'clusters', 3);
%         lyngby_km_plot_scat(centers, 'assign', assign, 'data', iris, 
%            'labels', kron([1 2 3],ones(1,50)))
%
%       See also LYNGBY, LYNGBY_KM_MAIN, LYNGBY_KM_PLOT_DIST.
%
% $Id: lyngby_km_plot_scat.m,v 1.2 2003/01/30 11:52:45 fnielsen Exp $


    % Sizes
    [K, dim] = size(Centers);

    % Check arguments
    if dim < 2
      error('Does not work with dimension < 2')
    end

    % Default properties
    assign = [];
    comp = [ 1 2 ];
    Data = [];
    labels = [];
    markerlist = '.ox+*sdv^<>ph';
    colorlist = 'mcrgbwk'; 
    
    % Parse properties
    n = 1;
    while n < length(varargin)
      arg = lower(varargin{n}); 

      if strcmp(arg, 'data') 
	n = n + 1;
	arg = varargin{n};
	if all(isreal(arg))
	  if size(arg,2) == dim
	    Data = arg;
	  else
	    error(['Argument to ''Data'' should have the same ' ...
		  'number of columns as ''Centers''.']);
	  end
	else
	  error(sprintf(['Argument with ''Data'' should be '...
		'reals.']));
	end
	 
      elseif strcmp(arg, 'assign') 
	n = n + 1;
	arg = varargin{n};
	if all(isreal(arg))
	  if max(arg) == K & min(arg) == 1
	    assign = arg(:);
	  else
	    error(['Argument to ''assign'' should be ' ...
		  'integers between 1 and K.']);
	  end
	else
	  error(sprintf(['Argument with ''Assign'' should be '...
		'integers.']));
	end

      elseif strcmp(arg, 'labels') 
	n = n + 1;
	arg = varargin{n};
	if all(isreal(arg))
	  if max(arg) == K & min(arg) == 1
	    labels = arg(:);
	  else
	    error(['Argument to ''Labels'' should be ' ...
		  'integers.']);
	  end
	else
	  error(sprintf(['Argument with ''Labels'' should be '...
		'integers.']));
	end

      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    if 0 | ~isempty(Data)
      [U, L, V] = svd(Data, 0);
    else
      [U, L, V] = svd(Centers, 0);
    end
    X = Centers * V(:,comp);
    
    cla
    % Plot every cluster center as a number
    for k = 1:K
      text(X(k, 1), X(k,2), int2str(k));
    end
    
    % Zoom out
    minX1 = min(X(:,1));
    maxX1 = max(X(:,1)); 
    minX2 = min(X(:,2));
    maxX2 = max(X(:,2));
    winX1 = maxX1 - minX1;
    winX2 = maxX2 - minX2;
    meanX1 = (minX1 + maxX1) / 2;
    meanX2 = (minX2 + maxX2) / 2;
    factor = 0.8;
    axis([ meanX1-winX1*factor  meanX1+winX1*factor ...
	  meanX2-winX2*factor meanX2+winX2*factor ]);

    if ~isempty(Data)
      Y = Data * V;
      hold on 
      if isempty(assign)
	plot(Y(:,comp(1)), Y(:,comp(2)), '.')
      else
	if isempty(labels)
	  for k = 1:K
	    i = find(assign == k);
	    c = colorlist(mod(k, length(colorlist))+1);
	    plot(Y(i,comp(1)), Y(i,comp(2)), [c '.']);
	  end
	else
	  for k = 1:K
	    iK = find(assign == k);
	    c = colorlist(mod(k, length(colorlist))+1);
	    for l = 1:K
	      iL = find(labels == l);
	      i = intersect(iL, iK);
	      m = markerlist(mod(l, length(markerlist))+1);
	      plot(Y(i,comp(1)), Y(i,comp(2)), [c m]);
	    end
	  end
	end
      end
    end

    for n = 1:2
      if comp(n) == 1
	label{n} = '1st principal component';
      elseif comp(n) == 2
	label{n} = '2nd principal component';
      elseif comp(n) == 3 
	label{n} = '3rd principal component';
      else 
	label{n} = sprintf('%dth principal component', n);
      end
    end
    xlabel(label{1})
    ylabel(label{2});
    title('Scattergram of cluster centers');













