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

% lyngby_km_error      - K-means, error 
%
%	function [e, mys, Sigmas, ps] = lyngby_km_error(X, a, ...
%           'PropertyName', 'PropertyValue')
%
%       Input:  X   Datamatrix, size: examples x variables. ie. if
%                   voxels are to be clustered the datamatrix should
%                   be voxels x time.
%               a   Assignment vector
%		Property:
%                  '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
%                  'TrainingSet'    Indices for the training set. The
%                                   first half of the dataset will be
%                                   used if it is not defined.
%                  'ValSet'         Indices for the validation set
%
% 	Output:	e   Error, negative log of p(x)
%
%       See also: lyngby_gmm_fit, lyngby_gmm_error

% cvs : $Id: lyngby_km_error.m,v 1.3 1998/05/11 14:46:03 fnielsen Exp $ 
%       $Revision: 1.3 $

    [rX, cX]  = size(X);

    % Default values
    clusterCenters = [];
    meanOrigin = 1;
    varianceType = 1;
    trainSet = [];
    valSet = [];
    reg = 0;
    
    n = 1;
    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 
      if strcmp(arg, 'meanorigin')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'meantrainingset')
	  meanOrigin = arg;
	elseif strcmp(arg, 'clustercenters')
	  meanOrigin = arg;
	else
	  error(sprintf('Type ''%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 = arg;
	elseif strcmp(arg, 'isotrop')
	  varianceType = arg;
	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

      elseif strcmp(arg, 'trainset')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if max(arg) <= size(X,1) & min(arg) > 0
	    trainSet  = arg;
	  else
	    error(sprintf(['The indices for ''TrainSet'' should ' ...
		'be larger than 0 and not larger than X. ' ...
		'Max index: %d, X: %d'], max(arg), size(X,1)));
	  end
	else
	  error(['The argument with ''TrainSet'' should be a ' ...
		'vector of indices.']); 
	end

      elseif strcmp(arg, 'valset')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if isempty(arg) | (max(arg) <= size(X,1) & min(arg) > 0)
	    valSet  = arg;
	  else
	    error(sprintf(['The indices for ''valSet'' should ' ...
		'be larger than 0 and not larger than X. ' ...
		'Max index: %d, X: %d'], max(arg), size(X,1)));
	  end
	else
	  error(['The argument with ''valSet'' should be a ' ...
		' vector containing indices.']); 
	end

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

    % Automatic training and validation set
    if isempty(trainSet)
      if isempty(valSet)
	index = floor(rX/2);
	trainSet = 1:index; 
	valSet   = (index+1):rX;
      else
	trainSet = lyngby_set_diff(1:rX, valSet);
      end
    elseif isempty(valSet)
      valSet = lyngby_set_diff(1:rX, trainSet);
    end
    if ~isempty(lyngby_set_and(trainSet, valSet))
      error('Training and test set should not be overlapping');
    end

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

    [mys, Sigmas, ps] = lyngby_gmm_fit(X(trainSet,:), ...
	'Assign', a(trainSet),  ...
	'MeanOrigin', meanOrigin, ...
	'ClusterCenters', clusterCenters, ...
	'VarianceType', varianceType, ...
	'Reg', reg);
    e = lyngby_gmm_error(X(valSet,:), mys, Sigmas, ps);











