function [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)

%SOM_PROBABILITY_GMM Probabilities based on a gaussian mixture model.
%
% [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)
% 
%   [K,P] = som_estimate_gmm(sM,D);
%   [pd,Pdm,pmd] = som_probability_gmm(D,sM,K,P);
%   som_show(sM,'color',pmd(:,1),'color',Pdm(:,1))  
%
%  Input and output arguments:
%   D    (matrix) size dlen x dim, the data for which the 
%        (struct) data struct,     probabilities are calculated
%   sM   (struct) map struct
%        (matrix) size munits x dim, the kernel centers
%   K    (matrix) size munits x dim, kernel width parameters
%                 computed by SOM_ESTIMATE_GMM
%   P    (matrix) size 1 x munits, a priori probabilities for each 
%                 kernel computed by SOM_ESTIMATE_GMM
%
%   pd   (vector) size dlen x 1, probability of each data vector in 
%                 terms of the whole gaussian mixture model
%   Pdm  (matrix) size munits x dlen, probability of each vector in 
%                 terms of each kernel
%   pmd  (matrix) size munits x dlen, probability of each vector to 
%                 have been generated by each kernel
%
% See also SOM_ESTIMATE_GMM.

% Contributed to SOM Toolbox vs2, February 2nd, 2000 by Esa Alhoniemi
% Copyright (c) by Esa Alhoniemi
% http://www.cis.hut.fi/projects/somtoolbox/

% ecco 180298 juuso 050100

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% input arguments
if isstruct(sM), M = sM.codebook; else M = sM; end
[c dim] = size(M);

if isstruct(D), D = D.data; end
dlen = size(D,1);

% reserve space for output variables
pd = zeros(dlen,1); 
if nargout>=2, Pdm = zeros(c,dlen); end
if nargout==3, pmd = zeros(c,dlen); end

% the parameters of each kernel
cCoeff = cell(c,1);
cCoinv = cell(c,1);
for m=1:c, 
  co = diag(K(m,:));
  cCoinv{m} = inv(co);
  cCoeff{m} = 1 / ((2*pi)^(dim/2)*det(co)^.5);
end

% go through the vectors one by one
for i=1:dlen, 

  x = D(i,:);
  
  % compute p(x|m)
  pxm = zeros(c,1); 
  for m = 1:c,
    dx     = M(m,:) - x;
    pxm(m) = cCoeff{m} * exp(-.5 * dx * cCoinv{m} * dx');
    %pxm(m) = normal(dx, zeros(1,dim), diag(K(m,:)));  
  end
  pxm(isnan(pxm(:))) = 0;
  
  % p(x|m)  
  if nargin>=2, Pdm(:,i) = pxm; end
  
  % P(x) = P(x|M) = sum( P(m) * p(x|m) )
  pd(i) = P*pxm; 
  
  % p(m|x) = p(x|m) * P(m) / P(x)
  if nargout==3, pmd(:,i) = (P' .* pxm) / pd(i); end
  
end


return; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% subfunction normal
%
% computes probability of x when mean and covariance matrix
% of a distribution are known

function result = normal(x, mu, co)

[l dim] = size(x);
coinv   = inv(co);
coeff   = 1 / ((2*pi)^(dim/2)*det(co)^.5);
diff   = x - mu;
result = coeff * exp(-.5 * diff * coinv * diff');