function Y = classify(X,mdl)
  for i = 1:length(mdl),
    res(i,:) = gaussian_prob(X, mdl(i).mu, mdl(i).Sigma);
    classidx(i) = mdl(i).class;
  end  
  [val,Y] = max(res);
  Y = classidx(Y);
  
function p = gaussian_prob(X, m, C);
  [D,N] = size(X);
  M = m'*ones(1,N);
  denom = 0.5*log(abs(det(C)));

  mahal = sum(((X-M)'* inv(C)).*(X-M)',2);
  p = -0.5*mahal - denom;