function [Y, c, Info] = lyngby_km_main(X, arg1, arg2, arg3, arg4, ...
    arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, ...
    arg15, arg16, arg17, arg18, arg19, arg20, arg21, arg22)

% lyngby_km_main       - Main function for K-means clustering 
%
%       function [Y, c, Info] = lyngby_km_main(X, ...
%           'PropertyName', 'PropertyValue')
%
%       Input:    X    Datamatrix, size: examples x variables. ie. if 
%                      voxels are to be clustered the datamatrix
%                      should be (voxels x time).
%
%       Property: Type            [ {median} | mean ]
%                 Standardization [ {None} | Std | Range ] Determines
%                                 the individual standardization
%                                 (normalization) of the variables
%                                 (the columns in the
%                                 datamatrix. 'Std' will standardize
%                                 with the standard deviation, 'Range'
%                                 with the difference max-min 
%                 Clusters        [ {10} | Integer ] Number of clusters.
%                 Init            [ {ReverseLog} | Linear |
%                                 UpperLinear | Random ] Initial
%                                 cluster centers determination. The
%                                 variables are sorted according to
%                                 max of xcorr or std of variables and
%                                 the initial centers are chosen from
%                                 this list. 
%                 DecayRate       Convergence control parameter, {0} <
%                                 DecayRate <= 1. Determines how the
%                                 clustering center converge. 
%                 Iterations      [ {20} | Integer ] Number of
%                                 iterations. 
%                 Variable        [ {time} | xcorr ] Clustering with
%                                 Cross correlation or time 
%                 Paradigm        Paradigm, the vector that is
%                                 used in the cross-correlation
%                                 with the datamatrix. This
%                                 variable needs to be defined if
%                                 the 'Variable' is 'xcorr'
%                 Components      [ {40} | integer ] Number of
%                                 cross-correlation components in the
%                                 analysis. Not used if 'Variabel' is
%                                 set to 'time'. Will max be set to
%                                 the number of columns in X 
%                 PositionWeight  Smoothing of the clustering. Weight
%                                 for the proximity part of the error
%                                 function. 
%
%       Output: Y       Cluster center matrix, size: Scans x 'Clusters' 
%               c       Assignment vector for the all voxels
%               Info    Shows the convergence of the Y's (array of
%                       lenght 'Iterations'). 
%
%       lyngby_km_main performs K-means or K-median clustering. The
%       number of clusters is specified with 'Clusters'. If 'Variable'
%       is 'time' then the datamatrix X will be used (directly) as the
%       input for the clustering algorithm. If 'Variable' is 'xcorr'
%       then the cross-correlation between the datamatrix and the
%       'Paradigm' will be used as input.
%
%       The individual variables (columns) in the datamatrix can be
%       scaled according to 'Standardization': With 'Std' the columns
%       are scaled to have equal standard deviation; with 'Range' the
%       difference between minimum and maximum in each column is used
%       to scale. Standardization should be used when the variables are
%       measured with different units or the interesting features
%       important for the discrimination lies in the variables with
%       low magnitude. When the centers are found they are scaled back
%       to the original space. 
%
%       'Init' determines how the cluster centers are initialized. For
%       all types of 'init' K specific objects (eg, voxels) are
%       selected (K corresponding to the number of clusters): For
%       'random' the initial cluster centers are initialized by
%       randomly picking K objects. For the other initialization
%       methods the selection is deterministic from sorted
%       objects. The sorting is either based on the standard deviation
%       of the original data or the maximum of the cross-correlation
%       function between the data and the paradigm. 'Linear' will
%       select with linear space though the sorted list of objects,
%       while 'reverselog' will select logarithmic through the list with
%       the most cluster centers picked from the objects with the
%       largest standard deviation or cross-correlation. 'UpperLinear'
%       will select from the top of the list.
%
%       Example: 
%         % K-means clustering of Fisher's iris data
%         load iris.txt
%         [Y,c,Info] = lyngby_km_main(iris,'type','mean','clusters',3);
%         figure,plot(Info),title('Convergence'),xlabel('Iteration');
%         C=zeros(3,3);for n=1:150,C(ceil(n/50.1),c(n))=C(ceil(n/50.1),c(n))+1;end
%         disp('Confusion matrix'), disp(C)
%
%       See also LYNGBY, LYNGBY_KM_CENTERSIM, LYNGBY_KM_PLOT_DIST,
%                LYNGBY_IKM_MAIN, LYNGBY_UI_KM_INIT, LYNGBY_XCORR.
%
% $Id: lyngby_km_main.m,v 1.28 2003/11/21 11:33:57 fnielsen Exp $ 


    % Default values
    decayRate  = 0;        
    Components = 40;
    init       = 1;
    iterations = 20;       
    K          = 10;       
    P          = [];
    positionWeight = 0;
    standard   = 1;
    Type       = 'median'; 
    Variable   = 'time';

    
    [M, N] = size(X);
    [rX, cX]  = size(X);
    
    n = 1;

    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 
      if strcmp(arg, 'type')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'mean')
	  Type = 'mean';
	elseif strcmp(arg, 'median')
	  Type = 'median';
	else
	  error(sprintf('Type ''%s'' not recognized', arg));
	end
	
      elseif strcmp(arg, 'standardization')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'none')
	  standard = 1;
	elseif strcmp(arg, 'std')
	  standard = 2;
	elseif strcmp(arg, 'range')
	  standard = 3;
	else
	  error(sprintf('Standardization ''%s'' not recognized', arg));
	end

      elseif strcmp(arg, 'clusters')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	K = arg;
	if K < 1 |  K > rX
	  error(sprintf(['The number of clusters should be ' ...
		'a value between 1 and the number of rows in X. ' ...
		'Clusters: %d; X: %dx%d'], K, size(X)));
	end
	
      elseif strcmp(arg, 'init')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if ~isstr(arg)
	  error('Argument to ''init'' should be a string');
	end
	if strcmp(arg, 'reverselog')
	  init = 1;
	elseif strcmp(arg, 'linear')
	  init = 2;
	elseif strcmp(arg, 'upperlinear')
	  init = 3;
	elseif strcmp(arg, 'random')
	  init = 4;
	else
	  error(sprintf('Argument with ''Init'' (''%s'') not recognized', arg));
	end

      elseif strcmp(arg, 'iterations')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	iterations = arg;
	if iterations < 1
	  error('The number of iterations should be a value larger than 1.')
	end
      
      elseif strcmp(arg, 'decayrate')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	decayRate = arg;
	if decayRate >= 1 | decayRate < 0
	  error('The DecayRate should be a value between 1 and 0 (0 < DecayRate <= 1).')
	end
	
      elseif strcmp(arg, 'positionweight')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	positionWeight = arg;
	if positionWeight <0
	  error('''PositionWeight'' should be a non-negative value.');
	end
	
      elseif strcmp(arg, 'variable')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'xcorr')
	  Variable = 'xcorr';
	elseif strcmp(arg, 'time')
	  Variable = 'time';
	elseif strcmp(arg, 'last results')
	  Variable = 'time';
	else
	  error(sprintf('Type ''%s'' not recognized', arg));
	end
	
      elseif strcmp(arg, 'paradigm')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	P = arg;
	if size(P,2) ~= size(X,2) & ~isempty(P)
	  error(sprintf(['The argument with ''Paradigm'' should ' ...
		'have the same number of columns as X. X: %dx%d, ' ...
		'Paradigm: %dx%d'], size(X), size(P)));
	end
	
      elseif strcmp(arg, 'components')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	Components = arg;
	if Components < 1
	  error('The number of components should be a value larger than 1.')
	end
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    if K > rX 
      error(sprintf(['The number of clusters should be less than ' ...
	    'the number of rows in X. Clusters: %d, X: %dx%d'], K, ...
	  size(X))); 
    end
    
    if strcmp(Variable, 'xcorr')
      if isempty(P)
	error(['P is empty: When ''Variable'' is ''xcorr'' ' ...
	    'then ''Paradigm'' should be defined']);
      end
    end

    if isempty(P)
      P = ones(rX,1);
    end
    
    if strcmp(Variable, 'xcorr')
      Components = min([Components size(X, 2)]);
    else
      Components = size(X, 2);
    end
    
    if positionWeight ~= 0
      
      dim  = lyngby_roi;
      xmin = dim(1,1);
      xmax = dim(1,2);
      ymin = dim(2,1);
      ymax = dim(2,2);
      zmin = dim(3,1);
      zmax = dim(3,2);
      
      SamplesX = (xmax-xmin+1);
      SamplesY = (ymax-ymin+1);
      SamplesZ = (zmax-zmin+1);
      
      if M ~= SamplesX*SamplesY*SamplesZ
	error([ 'This input volume is currently not allowed ' ...
	      '- must be square in 3 dims']);
      end

      index = 0 : SamplesX*SamplesY*SamplesZ-1;
      xpos  = lyngby_mod(index, SamplesX);
      rest  = round((index-xpos) / SamplesX);
      ypos  = lyngby_mod(rest, SamplesY);
      zpos  = round((rest-ypos) / SamplesY);
      
      R = sqrt(positionWeight) * [xpos' ypos' zpos'];
    end

      
    if strcmp(Variable,'xcorr')
      lyngby_log('Computing cross correlation');
      X = lyngby_xcorr(X', P', 'Components', Components, 'biased')';
    end

    % Initialize centers by K points in X with maximum variance
    if size(X,2) > 1
      if strcmp(Variable,'xcorr')
	[dummy, indx] = sort(max(X'));
      else
	[dummy, indx] = sort(std(X'));
      end
    else
      [dummy, indx] = sort(X');
    end

    % Standardization 
    if standard == 1
      % Do nothing
      standardValues = ones(1, size(X,2));
    elseif standard == 2
      lyngby_log('Scaling with standard deviation');
      standardValues = std(X);
      X = X ./ (ones(rX, 1) * standardValues + realmin);
    elseif standard == 3
      lyngby_log('Scaling with range');
      standardValues = max(X) - min(X);
      X = X ./ (ones(rX, 1) * standardValues + realmin); 
    end
    
    % Pad the additional columns for position dependent weighting 
    if positionWeight ~= 0
      X = [X R];
    end

    lyngby_log('Setting initial cluster positions');
    if init == 1
      lyngby_log('Initializing with reverse logspace');

      % There is a bug in the matlab version 5.5 of function
      % 'logspace' distributed with early versions of matlab 5.2. To
      % get a proper version of the logspace function the 5.2.1 patch
      % should be installed.
      
      % Reverse logspace
      % indx2 = rX+1 - round(logspace(log10(rX), 0, K));  
      indx2 = lyngby_set_unique(round(logspace(log10(rX), 0, K)));
      i = 1;
      while length(indx2) < K
	if indx2(i) ~= i
	  indx2 = [indx2(1:i-1) i indx2(i:end) ];
	end
	i = i + 1;
      end

    elseif init == 2
      lyngby_log('Initializing with linear space');
      indx2 = round(linspace(1, rX, K));                % Linear space 
    elseif init == 3
      lyngby_log('Initializing with upper linear space');
      step = 5;
      indx2 = [1 rX+1-(K-1)*step:step:rX];              % Upper Linear
    else
      lyngby_log('Initializing with random');
      indx2 = ceil(rX * rand(1, K));                    % Random
      indx2 = lyngby_set_unique(indx2);
      while length(indx2) < K
	indx2 = [indx2 ceil(rX * rand(1,K-length(indx2)))];
	indx2 = lyngby_set_unique(indx2);
      end
    end

    Y = X(indx(indx2),:);

    % First compute the matrix of input dot products
    if cX == 1
      xsum = (X.*X)';
    else
      xsum = sum((X.*X)');
    end
    D1 = xsum(ones(K,1), :)';

    Info = zeros(1, iterations);
    for t = 1:iterations


      % Distance matrix, MxK, 
      % Adding sum y^2 - x*y^iterations = sum(x-y)^2
      D = D1 + repmat(diag(Y*Y')',rX,1) - X*(2*Y)'; 
      % Profiling (the above is the slowest line in this function): 
      %   no speed enhancement: 
      %     ones(rX,1) * diag(Y*Y') -> repmat(diag(Y*Y')', rX, 1)
      %   50% speed enhancement:
      %     2*X*Y' -> X*(2*Y)'
      
      % c is the index matrix, ie. the assignment matrix
      if K == 1
	q = D';
	c = ones(1,rX);
      else
	[q, c] = min(D');
      end
      
      % Compute the centroid of points belonging to k'th center    
      da = 0;
      if strcmp(Type, 'mean')
	lyngby_log(sprintf('K-means iteration: %i / %i', t, iterations));
	
	% Mean
	for k = 1:K
	  a = Y(k,:);
	  objInCluster = sum(c==k);
	  if objInCluster > 1
	    Y(k, 1:Components) = Y(k, 1:Components) * decayRate + ...
		(1-decayRate) * mean(X(find(c==k),1:Components)); 
	    if positionWeight ~= 0
	      Y(k, Components+1:Components+3) = ...
		  Y(k, Components+1:Components+3) * decayRate ...
		  + (1-decayRate) * mean(X(find(c==k), ...
		  Components+1:Components+3)); 
	    end
	  elseif objInCluster == 1
	    Y(k, 1:Components) = Y(k, 1:Components) * decayRate + ...
		(1-decayRate) * X(find(c==k),1:Components);
	    if positionWeight ~= 0
	      Y(k, Components+1:Components+3) = ...
		  Y(k, Components+1:Components+3) * decayRate ...
		  + (1-decayRate) * X(find(c==k), ...
		  Components+1:Components+3); 
	    end
	  end
	  da = da + sum(abs(a-Y(k,:)));
	end
	
      else
	% Median
	lyngby_log(sprintf('K-mediod iteration: %i / %i', t, iterations));
	
	for k = 1:K
	  a = Y(k,:);
	  objInCluster = sum(c==k);
	  if objInCluster > 1
	    Y(k, 1:Components) = Y(k,1:Components) * decayRate + ...
		(1-decayRate) * median(X(find(c==k), 1:Components)); 
	    if positionWeight ~= 0
	      Y(k, Components+1:Components+3) = ...
		  Y(k, Components+1:Components+3) * decayRate ...
		  + (1-decayRate) * ...
		  median(X(find(c==k), Components+1:Components+3)); 
	    end
	  elseif objInCluster == 1
	    Y(k, 1:Components) = Y(k,1:Components) * decayRate + ...
		(1-decayRate) * X(find(c==k), 1:Components); 
	    if positionWeight ~= 0
	      Y(k, Components+1:Components+3) = ...
		  Y(k, Components+1:Components+3) * decayRate ...
		  + (1-decayRate) * X(find(c==k), ...
		  Components+1:Components+3); 
	    end
	  end
	  da = da + sum(abs(a-Y(k,:)));
	end
      end

      Info(t) = da; 
      if da == 0
	% If there is no change in the cluster centers we
	% should stop the iterations
	break;
      end
      
    end;

    Y = Y(:,1:Components) .* repmat(standardValues, size(Y,1), 1);
    c = c';












