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

% lyngby_gmm_error     - Gaussian mixture model, error
%
%	function e = lyngby_gmm_error(X, mys, Sigmas, ps)
%
%       Input:  X        Datamatrix, size: examples x
%                        variables. ie. if voxels are to be clustered
%                        the datamatrix should be voxels x time.
%               mys      Center of Gaussian mixture models
%               Sigmas   Covariance of Gaussian mixture models
%               ps       Mixing parameters
%
% 	Output:	e        Error, negative log of p(x)
%
%       Implements:
%           e = -sum(log(p(x)))
%
%           p(x) = sum_c( p_c |2pi Sigma_c|^(1/2) ...
%                               exp(-(x-my_c)' Sigma^(-1) (x-my_c)))

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

    [rX, cX]           = size(X);
    [rMys, cMys]       = size(mys);
    [rSigmas, cSigmas] = size(Sigmas);
    
    if rMys ~= length(ps) | rMys ~= size(Sigmas,1)
      error(sprintf(['mys, Sigmas and ps should have the same '...
	    'number of rows (length).\n' ...
	    'X: %dx%d, mys: %dx%d, Sigmas: %dx%d, ps: %dx%d'], ...
	  size(X), size(mys), size(Sigmas), size(ps)));
    end
    
    % Number of clusters
    K = rMys;
   
    P = zeros(K, rX);
    
    % Iterate of clusters
    for c = 1:K

      if ps(c) ~= 0
	% Centered datamatrix
	Xc  = X - ones(rX, 1) * mys(c,:);
      
	% Covariance
	Sigma = reshape(Sigmas(c,:), cX, cX);
	if isfinite(cond(Sigma))
	  isdSigma = 1 / sqrt(det(Sigma));
	  iSigma = inv(Sigma);
	  
	  if cX == 1
	    XSX = 0.5 * iSigma * (Xc .* Xc)';
	  else
	    XS  = Xc * iSigma;
	    XSX = 0.5 * sum((XS .* Xc)'); 
	  end
	  
	  P(c,:) = ps(c) * ((2*pi)^(-cX/2) * isdSigma) * exp(-XSX);
	end
      end
    end

    % Error as the negative log likelihood
    if size(P,1) == 1
      if any(P)
	e = - mean(log(P));
      else
	e = Inf;
      end
    else
      sumP = sum(P);
      if any(sumP)
	e = - mean(log(sumP));
      else
	e = Inf;  % log(0) produce warnings 
      end
    end
      














