% [mus, sigmas, weights] = trainGMM(X, mus, sigmas, weights, threshold, maxiterations, modeltype);
% 單撈磐 青紡 X惑俊輯 [mus, sigmas, weights]牢 頗扼固磐俊 狼竅咯 瀝秦柳 GMM闌  
% EM 舅絆府硫闌 撈儈茄 切嚼闌 肺弊快檔啊 THRESHOLD 藹焊促 利芭唱 
% MAXITERATIONS 館汗 寢薦啊 瞪 錠 鱉瘤 薦青.
% 去欽 己盒狼 傍盒魂俊 葷儈等 葛膽狼 屈怕綽 MODELTYPE俊輯 
% 膠納老等 竄困 青紡 趣籃 措阿 青紡肺 瘤瀝 啊瓷. 
% (肯傈 傍盒魂欄肺狼 犬厘且妨擱 sigmas啊 氦磐啊 酒聰扼 青紡肺 瀝秦廉具父 茄促.) 
%
% 澇仿 牢磊 :
% X綽 青撈 d瞞盔牢 漂隆 氦磐 單撈磐 青紡.
% mus : d*M 青紡肺 弊 青籃 雇膠媚 己盒俊 措茄 乞閉 氦磐.
% sigmas : d*M 青紡肺 弊 青籃 雇膠媚 去欽 己盒俊 措茄 盒魂 氦磐.
% weights : M-瞞盔 氦磐肺 去欽 啊吝摹.
% modeltype : 'scalarcov' 趣籃 'diagcov' 鞏磊凱.
% 免仿 牢磊 : 

function [mus, sigmas, weights] = trainGMM(X, mus, sigmas, weights, threshold, maxiterations, modeltype);
  
% Check input params and get dimensions:
  [d,numpoints] = size(X);
  M = length(weights);
  if (abs(sum(weights) - 1) < 0.00001) == 0
     error('trainGMM: Weights should sum to 1!');
  end
  
  if (size(mus) == [d,M] & size(sigmas) == [d,M]) == 0
     error('trainGMM: Weights should sum to 1!');
  end
 
  if strcmp(modeltype,'scalarcov')
    xvars = zeros(1,numpoints);
  elseif strcmp(modeltype,'diagcov')
    xvarmat = zeros(d,numpoints);
  else
    error('Unknown model type!');
  end

% Plot data:
  % hold off
  % plot(X(1,:),X(2,:),'r.')
  % hold on

% Plot initial component centroids:
  %  plot(mus(1,:),mus(2,:),'bo')

% Calculate initial likelihoods:
  [loglikevec, pmat] = loglikeGMM(X, mus, sigmas, weights);
  loglike = sum(loglikevec);
  disp('Log likelihood at initialisation: '), disp(loglike);

% Iterate until convergence:
  for iteration = 1:maxiterations
    disp('Iteration:'),disp(iteration);

% E step: calculate posteriors under current parameter values:
    for n = 1:numpoints
      pmat(:,n) = pmat(:,n)/sum(pmat(:,n));
    end
    disp(pmat(:,1:5))

% M step: update parameter values to maximise posteriors:
    for m = 1:M;
% Mean update:
      psum = sum(pmat(m,:));
      mus(:,m) = (X * pmat(m,:)')/psum;
% Variance update:
      if strcmp(modeltype,'scalarcov')
	for n = 1:numpoints
          meansub = X(:,n)-mus(:,m);
          xvars(n) = meansub'*meansub;
	end
	sigmas(:,m) = ones(d,1)*(xvars * pmat(m,:)')/(d*psum);
      elseif strcmp(modeltype,'diagcov')
	for n = 1:numpoints
          meansub = X(:,n)-mus(:,m);
          xvarmat(:,n) = meansub.*meansub;
	end
	sigmas(:,m) = xvarmat * pmat(m,:)' / psum;       
      end
% Weight update:
      weights(m) = psum/numpoints;
    end

% Plot new component centroids:
  drawnow
  plot(mus(1,:),mus(2,:),'bo') % 頗扼固磐 彌利拳 苞瀝闌 焊咯林扁 困秦輯
  disp(weights)
  disp(mus)
  disp(sigmas)

% Calculate new log likelihood:
    [newloglikevec, pmat] = loglikeGMM(X, mus, sigmas, weights);
    newloglike = sum(newloglikevec);
    disp('Log likelihood: '), disp(newloglike);

% Test for convergence:
    if (newloglike < loglike)
      error('Log likelihood decreased!')
    end
    if ((newloglike - loglike)/numpoints < threshold)
      disp('Converged.')
      break;
    else 
      loglike = newloglike;
    end;

  end
  disp(' '),disp(' '),disp(' ')