function [memberships,logp] = gaussian_mixture_estep(data,params);
%  function [memberships, logp] = gaussian_mixture_estep(data,params);
%
%  Estimates membership probabilities (E-step) for a Gaussian mixture model 
%  with d variables, given a model specified by the data structure "params".  
%
%  INPUTS
%  data: N x d real-valued data matrix
%  params:  2d array of structures of the form: 
%           params(k).weight = weight of component k
%           params(k).mean = d-dimensional mean vector for kth component 
%           params(k).covariance = d x d covariance vector for kth
%           component
% 
%  OUTPUTS
%  memberships = n x K matrix of component memberships
%                (note that memberships for each row sum to 1)
%
%
%                                        Padhraic Smyth, CS 274A, Winter 2009


% Ideally need to add some error-checking here....

[n d] = size(data);
K = size(params,2);
 

% precompute log-weights and log-conditionals to save time
for k=1:K  %  for each mixture component....
    params(k).logweight = log(params(k).weight); 
end

% compute the n x K matrix of log conditional probabilities, log p(x_ij|k):
for k=1:K 
    logp_xi_given_k(:,k) = log(pgauss_chol(data,params(k).mean, params(k).covariance));
end 

 
% compute the matrix logp(xi,k), an n x K matrix of joint probabilities
% (sum over the columns (variables) of logp_xi_given_k.....)
for k=1:K
    logp_xi_and_k(:,k) =  logp_xi_given_k(:,k)  + repmat(params(k).logweight,n,1);
end
 
% Identify and remove the max element of each row to prevent underflow in exponentiation  
    c = max(logp_xi_and_k')';
% compute a matrix of numbers proportional (within c) to the joint probabilities p(xi,k) 
   p_xi_and_k = exp(logp_xi_and_k - repmat(c,1,K));   
   %   p_xi_and_k = exp(logp_xi_and_k);  
   
% sum out over the K components (sum each row) to get numbers proportional (within c) to p(x_i)
    p_xi = sum(p_xi_and_k')';  
% compute the relative component probabilities p(k|xi), (the c factor cancels out top and bottom)
    memberships = p_xi_and_k./repmat(p_xi,1,K);

% Compute the log-likelihood, making sure to add back in the factor of c for
% each row that was removed earlier 
    logp = sum(log(p_xi)+c);