function [Centers, Assign, Within, Property] = lyngby_ikm_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_ikm_main      - Iterative K-means clustering  
%
%       function [Center, Assign, Within, Prop] = lyngby_ikm_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         List of clusters {1:5}
%                 Runs             [ {1} | Integer ] Number of runs 
%                 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:   Center     Cluster center matrix, size: cluster x variables
%                 Assign     Assignment vector for the all objects
%                 Within     Trace of Within matrix cluster x runs  
%                 Property   Properties
%
%       "Iterative K-means" iterate over a list of number of clusters
%       each time calling lyngby_km_main. The intra-cluster variance
%       is calculate for each configuration. 
%
%       See also LYNGBY, LYNGBY_KM_MAIN, LYNGBY_KM_WITHIN,
%                LYNGBY_IKM_PLOT_W, LYNGBY_IKM_BIC.   
%
% $Id: lyngby_ikm_main.m,v 1.5 2003/02/05 18:53:27 fnielsen Exp $ 


    % Default values
    Clusters   = 1:5;
    Components = 40;
    decayRate  = 0;        
    init       = 'random';
    iterations = 20;       
    positionWeight = 0;
    P          = [];
    Runs       = 1;
    standard   = 'none';
    Type       = 'mean'; 
    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 = 'none';
	elseif strcmp(arg, 'std')
	  standard = 'std';
	elseif strcmp(arg, 'range')
	  standard = 'range';
	else
	  error(sprintf('Standardization ''%s'' not recognized', arg));
	end

      elseif strcmp(arg, 'clusters')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	Clusters = arg;
	if any(Clusters < 1) |  any(Clusters > rX)
	  error(sprintf(['The number of clusters should be ' ...
		'a list with values between 1 and the number of rows in X. ' ...
		'Clusters: %d; X: %dx%d'], Clusters, size(X)));
	end
	
      elseif strcmp(arg, 'runs')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	Runs = arg;
	if Runs < 1
	  error(sprintf(['Runs should be ' ...
		'a value larger than 1. ' ...
		'Runs: %d; X: %dx%d'], Runs, 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 = 'reverselog';
	elseif strcmp(arg, 'linear')
	  init = 'linear';
	elseif strcmp(arg, 'upperlinear')
	  init = 'upperlinear';
	elseif strcmp(arg, 'random')
	  init = 'random';
	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)
	  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(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

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

    if strcmp(Variable, 'xcorr')
      Components = min([Components size(X, 2)]);
    else
      Components = size(X, 2);
    end
    
    Within = Inf * ones(length(Clusters), Runs);
    Centers = [];
    
    for n = 1:length(Clusters)
      for r = 1:Runs
	lyngby_log(sprintf('Starting run %i using %i clusters',r,Clusters(n)));
	[aCenters, aAssign, aInfo] = ...
	    lyngby_km_main(X, ...
	    'Type', Type, ...
	    'Standardization', standard, ...
	    'Clusters', Clusters(n), ...
	    'Init', init, ...
	    'DecayRate', decayRate, ...
	    'PositionWeight', positionWeight, ...
	    'Iterations', iterations, ...
	    'Variable', Variable, ...
	    'Paradigm', P, ...
	    'Components', Components);
	Within(n,r) = lyngby_km_within(X, aCenters, aAssign, ...
	    'type', 'trace'); 
	if r == 1 | min(Within(n,1:r-1)) > Within(n,r)
	  if n == 1
	    indices = 1:Clusters(1);
	  else
	    start = sum(Clusters(1:n-1));
	    indices = start+1:start+Clusters(n);
	  end
	  Centers(indices,:) = aCenters;
	  Assign(:,n) = aAssign;
	end
      end
    end

    Property.Type = Type;
    Property.Standardization = standard;
    Property.Clusters = Clusters;
    Property.Runs = Runs;
    Property.Init = init;
    Property.DecayRate = decayRate;
    Property.Iterations = iterations;
    Property.Variable = Variable;
    Property.Paradigm = P;
    Property.Components = Components;
    Property.PositionWeight = positionWeight;
    




