function p = pgauss_chol(x, m, s)
%pgauss_chol	gaussian probability via cholesky factorization
% 
% p = pgauss_chol(x, m, s)
% * Compute the probability of each row vector in x assuming
% data is from a multivariate normal distribution with mean m and
% covariance matrix s. The mean may be a row or a column vector.
% * This is a recasting of Mike Burl's original code for this purpose;
% that code does not perform well for large nprob.
% * Emphasize that if xlen=1, s is the variance, not standard deviation.
% * The bulk of the computation is in finding, essentially,
%     (x-m) s^(-1) (x-m)' 
%   = (x-m) (r' r)^(-1) (x-m)' 
%   = (x-m) r^(-1) r^(-1)' (x-m)' = || (x-m) r^(-1) ||_2
% where r is the Cholesky factor of s, and ||_2 is the 2-norm of the
% rows of the enclosed matrix.  The original version found the entire
% product which means forming an nprobXnprob matrix, unworkable if
% nprob is large.  This approach never finds a matrix of size larger
% than that of x.
% 
% Inputs:
%   real x[nprob,xlen];
%   real m[xlen];
%   real s[xlen,xlen];  -- positive-definite
% 
% Outputs:
%   real p[nprob,1];
% 
% See Also:  pgauss

% 
% Error checking
% 
if all(nargin  ~= [3]), error ('Bad input arg number'); end
% if all(nargout ~= [0 1]), error ('Bad output arg number'); end  
% check sizes
[nprob, xlen] = size(x); % number of probs to find, length of one sample
if nprob == 0, p = []; return; end;
if (length(m) ~= xlen), error('x and m not conformant'); end;
if any(size(s) ~= [xlen xlen]), error('m and s not conformant'); end;

%
% Computation
% 
m = m(:)'; % force a row
r = chol(s); % square root of s, upper-triangular

% the log of the scale factor, identical for all samples
% NB prod(diag(r)) = sqrt(det(s))
l_factor = -(xlen/2) * log(2*pi) - sum(log(diag(r)));

% special case because sum(vector) will sum the whole thing to a scalar
if xlen > 1,
   p = exp(l_factor - 0.5 * sum(((x-m(ones(nprob,1),:)) * inv(r))'.^2))';
else
   p = exp(l_factor - 0.5 * ((x-m) / r).^2 );
end;
return;