function [mask P M S] = mvgmm_classifier (x,K)

%PUBLISHED UNDER THE CREATIVE COMMONS LICENSE (SEE MY WEBPAGE FOR MORE DETAILS)

%INPUT:
%x is the input vector where clustering is performed, x in R^d
%K is the number of clusters

%OUTPUT:
%mask is the classification of points of x => mask(i) is the category of point x(i), categories in {1,2,3,...,K}
%P: priors of mixture model
%M: centers of clusters
%S: covariances of clusters

%to provide clusters' centers' initialization set the global variable MU_INIT
%MU_INIT is an Kxd matrix
global MU_INIT;

EPSILON = 1E-06;
ITER    = 20;    % EM LOOP

N = size (x,1); d = size (x,2);

P0 = zeros (K,1);
M0 = zeros (K,d);
S0 = zeros (d,d,K);

ti = zeros (K,N);

Emax = inf; Enew = inf;

if isempty (MU_INIT)
    %INITIALIZE WITH K-MEANS
    mask = kmeans (x, K, 'emptyaction','singleton','display','off',...
        'start','cluster','replicates',20);

    %INITIALIZE PARAMETERS ACCORDING TO MASK
    for i=1:K
        ind = find (mask == i);
        if isempty (ind)
            P0 (i)     = EPSILON;
            S0 (:,:,i) = 1000*eye(d);
        else
            P0 (i)     = length(ind) / N;
            M0 (i,:)   = mean (x(ind,:));
            S0 (:,:,i) = cov(x(ind,:));
            %IF COVARIANCE IS NOT POSITIVE DEFINITE (=NEGATIVE OR ZERO
            %EIGENVALUES) => CLUSTER IS USELESS => LARGE COVARIANCE
            [veig e] = eig (S0(:,:,i));
            if ~isempty (find(diag(e) < 0))
                P0 (i)    = EPSILON;
                S0(:,:,i) = 1000 * eye(d);
            else
                e = diag(e); ind=find(e<1.); e(ind)=1.; e=diag(e);
                S0(:,:,i) = veig*e*veig';
            end
        end
    end
else
    M0 = MU_INIT;
    for i=1:K
        S0(:,:,i) = eye (d);
    end
    P0 = (1/K)*ones(K,1);
end

%-

%EM LOOP
Eold = inf;
for em_loop=1:ITER
    %E-STEP
    px = zeros (N,1);
    for i=1:K
        px = px + P0(i)*mvnpdf (x,M0(i,:),S0(:,:,i));
    end
    px(find(px == 0)) = EPSILON;
    for i=1:K
        S1 = inv (S0(:,:,i));
        %POSTERIOR
        ti(i,:) = P0(i)*mvnpdf(x,M0(i,:),S0(:,:,i))./px;

        %CHECK FOR ZERO VALUES
        ti(i,find(ti(i,:)<EPSILON)) = EPSILON;
    end

    %NEGATIVE LOGARITHMIC LIKELIHOOD (MINIMIZATION)
    px = sum(ti,2); Eold = Enew; Enew = - sum(log(px));
    %CHECK FOR CONVERGENCE
    if Eold - Enew < EPSILON
        break;
    end
    
    P = P0;
    M = M0;
    S = S0;

    %M-STEP
    for i=1:K
        P0(i) = sum(ti(i,:))/N;

        M0(i,:) = sum(repmat(ti(i,:)',1,d).*x)/sum(ti(i,:));

        z = (x-repmat(M0(i,:),N,1))'; R = zeros(size(S0(:,:,i)));
        for j=1:N
            R = R + ti(i,j)*z(:,j)*z(:,j)';
        end
        S0(:,:,i) = R/sum(ti(i,:));
        [veig e] = eig (S0(:,:,i)); e = diag(e);
        if ~isempty (find(diag(e) < 0))
            P0 (i)    = EPSILON;
            S0(:,:,i) = 1000 * eye(d);
        end
    end    
end

%NORMALIZATION OF PRIORS
P = P / sum(P);

%FIND MASK
px = zeros (N,1);
for i=1:K
    px = px + P(i)*mvnpdf (x,M(i,:),S(:,:,i));
end
px(find(px == 0)) = EPSILON;
for i=1:K
    %POSTERIOR
    ti(i,:) = P(i)*mvnpdf(x,M(i,:),S(:,:,i))./px;

    %CHECK FOR ZERO VALUES
    ti(i,find(ti(i,:)<EPSILON)) = EPSILON;
end
[q mask] = max (ti,[],1);