function [params,memberships] =  gaussian_mixture(data,K,initialization_method,RSEED,niterations,plotgap,initparams)
% function [params,memberships] =  gaussian_mixture(data,K,initialization_method,RSEED,niterations)
%
% fits a Gaussian mixture model via EM to d-dimensional data using K
% mixture components. The data are contained in the N x d "data" matrix. 
%  
%
% INPUTS
%  data: N x d real-valued data matrix
%  K: number of clusters (mixture components)
%  initialization_method: 1 for kmeans, 2 for random memberships
%
% OUTPUTS
%  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
%
%  memberships: n x K matrix of probability memberships for "data"
%
% Note 1: Interpretation of Params and Memberships:
%    - params(k).weight is the probability that a randomly selected row
%         belongs to component (or cluster) i (so it is "cluster size")
%    - memberships(i,k) = p(cluster k | x) which is the probability
%         (computed via Bayes rule) that vector x was generated by cluster
%         k, according to the "generative" probabilistic model.
%
%
%                                        Padhraic Smyth, CS 274A, Winter 2009
 
 
[n d] = size(data);
EPS = 0.00001;     % if the relative change in log-likelihood from
                   % one iteration to the next is less than EPS, then
                   % halt the iterations and declare convergence.
MIN_ITERATIONS = 3;  % ensures some minimum number of iterations are run
MAX_ITERATIONS = niterations;  % the maximum number of iterations to perform
                               % (useful to prevent "runaway" slow-converging
                               % runs, particularly on large data sets). 
plotflag=1;
if nargin<6
    plotgap=1;
end

% Ideally should add some checking here that EPS>0, that
% MIN_ITERATIONS < MAX_ITERATIONS, that the input arguments have
% valid dimensionalities (e.g., that n>1, d>1, etc).

% select initial parameter values for the mixture model
if nargin<7
    %initialization_method = 1;  % 1 = kmeans, 2 = random initialization 
    fprintf('Running initialization....\n\n');
    memberships = gaussian_mixture_initialize(data,K,RSEED,initialization_method);
else
    params =  initparams;
end


 
% begin EM iterations until convergence
iteration = 0;
threshold_condition=1; 


while (threshold_condition == 1)
      iteration = iteration + 1;
       
        
      % Mstep: calculate the new parameters, given the memberships from the
      % E-step
      [newparams] = gaussian_mixture_mstep(data,memberships);
      
      % Estep: calculate component memberships     
      [memberships,Lnew] = gaussian_mixture_estep(data,newparams);
      fprintf('Iteration %d: Log-likelihood per data point = %10.4f \n',iteration,Lnew/n);    
      % Note that the log-likelihood Lnew cannot decrease in an EM
      % algorithm, so Lnew/n should never decrease.
      if plotflag==1 
          if iteration==1
               plotdata(data);
               pause;
               plot_gaussians(data, newparams,1,2,Lnew/n,iteration,memberships); 
               pause;
          else
              if mod(iteration,plotgap)==0
              plot_gaussians(data, newparams,1,2,Lnew/n,iteration,memberships); 
              end
          end     
      end 
      
      % store the change in likelihood between iteration 1 and iteration 2
      % to provide a scale to determine convergence  
      if (iteration == 2)
         Ldelta2 = Lnew - Lold;
      end
      
      % check for convergence (within a specified range of minimum and
      % maximum number of iterations)
		if iteration>MIN_ITERATIONS
         Ldelta = Lnew - Lold;
         if Ldelta<(Ldelta2*EPS) | iteration>(MAX_ITERATIONS-1)
         	threshold_condition = 0;
         end
        end
      
      % update...
      Lold = Lnew;
      params = newparams;  
       logl(iteration) = Lnew;
end 
   

if plotflag==1
    figure;
    plot(1:iteration,logl./n,'r-');
    title('Log-Likelihood by Iteration');
    xlabel('Iteration Number');
end