function S = lyngby_km_centersim(X, Centers, assign, varargin) 

% lyngby_km_centersim  - K-means center similarity
%
%	function S = lyngby_km_centersim(X, Centers, assign)
%
%       Input:    X          Datamatrix, size: objects x variables 
%                 Centers    Cluster centers, size: clusters x
%                            variables 
%                 assign     Assignment vector, size: objects x 1 
%
%       Property: Similarity [ {InverseEuclidean} | corrcoef ]
%
%       Output:   S          Similarity matrix, size: objects x 1 
%
%       This function will calculate the similarity between the
%       cluster centers and their assigned data points from the X
%       matrix. 
%       
%       If 'similarity' is 'InverseEuclidean' the similarity is
%       calculated as the inverse of the Euclidean distance:
%           1/sum( (X-center).^2 )
%       For 'CorrCoef' as 'similarity' the similarity is computed as
%       the correlation coefficient, that is between -1 and +1 as
%       implemented in the corrcoef Matlab function
%
%       See also LYNGBY, LYNGBY_KM_MAIN, LYNGBY_KM_PLOT_DIST.
%
% $Id: lyngby_km_centersim.m,v 1.3 2003/01/30 12:03:20 fnielsen Exp $


    % Sizes
    [rX, cX] = size(X);
    [rAssign, cAssign] = size(assign);
    [rCenters, cCenters] = size(Centers);

    nClusters = rCenters;

    % Check arguments
    if nargin < 3
      error('Too few input arguments');
    end
    if rX ~= rAssign
      error('X and assign should have the same number of rows');
    end
    if cX ~= cCenters
      error('X and Centers should have the same number of columns');
    end

    
    % Default values
    similarity = 'inverseeuclidean'; 


    % Parse properties
    n = 1;
    while n < length(varargin)
      arg = lower(varargin{n});

      if strcmp(arg, 'similarity') 
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg)
	  if strcmp(arg, 'inverseeuclidean') | strcmp(arg, 'corrcoef')
	    similarity = arg;
	  else
	    error(['Argument to ''Similarity'' should be ' ...
		  ' ''InverseEuclidean'' or ''CorrCoef''.']);
	  end
	else
	  error(sprintf(['Argument with ''Similarity'' should be '...
		'a string']));
	end

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

    
    
    S = zeros(rX, 1);
    
    for n = 1:nClusters
      index = find(assign == n); 
      if strcmp(similarity, 'inverseeuclidean')
	dist = sum((X(index,:) - Centers(n*ones(1,length(index)),:)).^2,2);
	indexZero = find(~dist);
	if ~isempty(indexZero)
	  indexNonzero = setdiff(1:length(dist), indexZero);
	  S(index(indexZero),1) = Inf;
	  S(index(indexNonzero),1) = 1./dist(indexNonzero);
	else
	  S(index,1) = 1./dist;
	end
      elseif strcmp(similarity, 'corrcoef')
	Xm = X(index,:) - repmat(mean(X(index,:),2), 1, size(X,2));
	Cm = Centers(n,:) - repmat(mean(Centers(n,:)), 1, size(X, 2));
	Cm = repmat(Cm, length(index), 1);
	S(index, 1) = sum(Xm .* Cm, 2) ./ (sqrt(sum(Xm.*Xm,2) .* sum(Cm.*Cm,2)));
      else
	error('Internal error')
      end
    end

