function [mys, Sigmas, ps] = lyngby_gmm_fit(X, arg1, arg2, arg3, ...
    arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, ...
    arg14, arg15, arg16, arg17, arg18, arg19, arg20, arg21, arg22)

% lyngby_gmm_fit       - Gaussian mixture model, fitting 
%
%	function [mys, Sigmas, ps] = lyngby_gmm_fit(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:
%                  'Assign'         Assignment vector, predefine
%                                   cluster membership
%                  'MeanOrigin'     [ {MeanTrainingSet} | ClusterCenters ]
%                                   Specifices who the center of the
%                                   gaussians should be found.
%                  'ClusterCenters' Centers for the clusters. This
%                                   should be defined if 'MeanOrigin'
%                                   is set to 'ClusterCenters'
%                  'VarianceType'   [ {Anisotrop} | Isotrop ]
%                  'Reg'            {0} Regularization
%
% 	Output:	mys      Centers of Gaussian clusters
%               Sigmas   Covariance of Gaussian clusters
%               ps       Mixing parameters

% cvs : $Id: lyngby_gmm_fit.m,v 1.1 1998/05/11 14:43:24 fnielsen Exp $ 
%       $Revision: 1.1 $

    [rX, cX]  = size(X);

    % Default values
    clusterCenters = [];
    meanOrigin = 1;
    varianceType = 1;
    reg = 0;
    assign = [];
    
    n = 1;
    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 
      if strcmp(arg, 'assign')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~isstr(arg)
	  if length(arg) == rX
	    assign = arg;
	  else
	    error(sprintf(['Argument with ''assign'' should be '...
		  'of the same size as the datamatrix X.\n' ...
		  'X: %dx%d, assign: %dx%d'], size(X), size(arg)));
	  end
	else
	  error('Argument to ''assign'' should be a vector');
	end

      elseif strcmp(arg, 'meanorigin')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'meantrainingset')
	  meanOrigin = 1;
	elseif strcmp(arg, 'clustercenters')
	  meanOrigin = 2;
	else
	  error(sprintf('MeanOrigin ''%s'' not recognized', arg));
	end
	
      elseif strcmp(arg, 'clustercenters')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if size(arg,2) == size(X,2) 
	    clusterCenters  = arg;
	  else
	    error(sprintf(['ClusterCenters'' should have the ' ...
		  'same number of columns as X. ' ...
		  'Cluster center: %dx%d, X: %dx%d'], ...
		size(arg), size(X)));
	  end
	else
	  error(['The argument with ''TrainSet'' should be a ' ...
		'vector of indices.']); 
	end

      elseif strcmp(arg, 'variancetype')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'anisotrop')
	  varianceType = 1;
	elseif strcmp(arg, 'isotrop')
	  varianceType = 2;
	else
	  error(sprintf('VarianceType ''%s'' not recognized', arg));
	end

      elseif strcmp(arg, 'reg')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~isstr(arg)
	  if arg >= 0
	    reg = arg;
	  else
	    error('Argument to ''Reg'' should be a larger than 0');
	  end
	else
	  error('Argument to ''Reg'' should be a float');
	end

      else
	error('Invalid property');
      end
      n = n + 1;
    end

    if meanOrigin == 2
      if isempty(clusterCenters)
	error(['When the mean of the gaussian is taken to be ' ...
	      'the cluster centers ''ClusterCenters'' should be ' ...
	      'defined']);
      end
    end

    % Number of clusters
    K = max(assign);
    
    mys = zeros(K, cX);
    Sigmas = zeros(K, cX*cX);
    ps = zeros(K, 1);
    
    % Iterate of clusters
    for c = 1:K
      ps(c) = sum(assign==c) / length(assign);

      indices = find(assign==c);
      if ~isempty(indices) & ~(length(indices) == 1 & reg == 0)
	% Mean
	if meanOrigin == 1
	  if length(indices) == 1
	    my = X(indices,:);
	  else
	    my = mean(X(indices,:));
	  end
	elseif meanOrigin == 2
	  my = clusterCenters(c,:);
	end
	
	% Variance
	if varianceType == 1
	  % Anisotrop variance	  
	  if length(indices) ~= 1
	    Sigma = cov(X(indices,:)) + reg * eye(cX);
	  else
	    Sigma = reg * eye(cX);
	  end
	elseif varianceType == 2	  
	  % Isotrop variance
	  if length(indices) ~= 1
	    Sigma = trace(cov(X(indices,:))) / cX + reg;
	  else
	    Sigma = reg;
	  end
	end

	mys(c,:) = my;
	if varianceType == 1
	  Sigmas(c,:) = Sigma(:)';
	else
	  Sigmas(c,:) = reshape(Sigma * eye(cX), 1, cX*cX);
	end
      end
      
    end














