function out = gaussianLog(data, gaussianParam, gConst);
% gaussianLog: Multi-dimensional log Gaussian propability density function
%	Usage: out = gaussianLog(data, gaussianParam, gConst)
%	data: d x n data matrix, representing n data vector of dimension d
%	mu: d x 1 vector
%	sigma: d x d matrix
%		d x 1 vector
%		1 x 1 scalar
%	out: 1 x n vector
%
%	See also "gaussian".

%	Roger Jang, 20050327

if nargin<1, selfdemo; return; end
[dim, dataNum]=size(data); 

dataMinusMu = data-repmat(gaussianParam.mu, 1, dataNum);
if prod(size(gaussianParam.sigma))==1	% identity covariance matrix times a constant for each Gaussian
	if nargin<3
		out = -sum(dataMinusMu.*dataMinusMu/gaussianParam.sigma, 1)/2-dim/2*log(2*pi)-dim/2*log(gaussianParam.sigma);
	else
		out = -sum(dataMinusMu.*dataMinusMu/gaussianParam.sigma, 1)/2+gConst;
	end
elseif prod(size(gaussianParam.sigma))==dim	% diagonal covariance matrix for each Gaussian
	if nargin<3
		out = -sum(dataMinusMu./repmat(gaussianParam.sigma, 1, dataNum).*dataMinusMu, 1)/2-dim/2*log(2*pi)-log(prod(gaussianParam.sigma))/2;
	else
		out = -sum(dataMinusMu./repmat(gaussianParam.sigma, 1, dataNum).*dataMinusMu, 1)/2+gConst;
	end
else	% full covariance matrix for each Gaussian
	invCov = inv(gaussianParam.sigma);			% For repeated invocation of this function, this step should be moved out of this function
	if nargin<3
		out = -sum(dataMinusMu.*(invCov*dataMinusMu), 1)/2-dim/2*log(2*pi)-log(det(gaussianParam.sigma))/2;
	else
		out = -sum(dataMinusMu.*(invCov*dataMinusMu), 1)/2+gConst;
	end
end

% ====== Self demo
function selfdemo
x = linspace(-10, 10);
gaussianParam.mu = 0;
gaussianParam.sigma = 0.1;

y1 = log(gaussian(x, gaussianParam));
y2 = gaussianLog(x, gaussianParam);
diff = abs(y1-y2);

subplot(2,1,1); plot(x, y1, x, y2);
subplot(2,1,2); plot(x, diff);