function [mask, P, M, S, V] = mvtmm_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
%V: degrees of freedom 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);
V0 = 5 + randn(K,1);

ti = zeros (K,N); u = zeros (size(ti)); dM = zeros (N,1);

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);
            V0(i)      = 20;
        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)); e = diag(e);
            if ~isempty (find(e < 0))
                P0 (i)    = EPSILON;
                S0 (:,:,i) = 100 * eye(d);
            else
                ind1=find(e<1.); e(ind1)=1.;  %ELIMINATE SMALL EIGENVALUES
                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
for em_loop=1:ITER
    %E-STEP
    px = zeros (N,1);
    for i=1:K
        px = px + P0(i)*mvtpdf (x,M0(i,:),S0(:,:,i),V0(i));
    end
    px(find(px == 0)) = EPSILON;
    for i=1:K
        %POSTERIOR
        ti(i,:) = P0(i)*mvtpdf(x,M0(i,:),S0(:,:,i),V0(i))./px;

        %U WEIGHT (SEE STUDENT t DISTRIBUTION DOCUMENTATION FOR U)
        S1 = inv (S0(:,:,i));
        for j=1:N
            dM(j) = (x(j,:) - M0(i,:)) * S1 * (x(j,:) - M0(i,:))';
        end
        u(i,:) = (V0(i) + d) ./ (V0(i) + dM);

        %CHECK FOR ZERO VALUES
        u (i,find(u(i,:)<EPSILON))  = EPSILON;
        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;
    V = V0;
    
    %M-STEP
    for i=1:K
        P0(i) = sum(ti(i,:))/N;

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

        z = (x-repmat(M0(i,:),N,1))'; R = zeros(size(S0(:,:,i)));
        for j=1:N
            R = R + ti(i,j).*u(i,j)*z(:,j)*z(:,j)';
        end
        S0(:,:,i) = R/sum(ti(i,:));
        %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);
            V0(i)     = 10;
        else
            V0(i) = get_dof (V0(i),d,u(i,:),ti(i,:));
        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)*mvtpdf(x,M(i,:),S(:,:,i),V(i))./px;

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

%
% -
%

%
%function:get_dof
%
function [new_v] = get_dof (v,p,u,t)

v = (v+p)/2; n = sum(t,2);
Q = 1+sum(t.*(log(u)-u))/n+psi(v)-log(v);

fun_z = @(x)((-psi(x/2)+log(x/2)+Q)^2);

[new_v val] = fminbnd (fun_z, 1E-05, 1000);
if new_v<1E-05 new_v=1E-05; end

%
%function:mvtpdf
%
function [p] = mvtpdf (x,m,s,v)

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

p = zeros (N,1);

Dt = det(s)^.5; S = inv(s);
for i=1:N
    t       = (x(i,:)-m)';
    d_mahal = t' * S * t;
    d_mahal = (1 + d_mahal / v)^((v+d)/2);
    p(i)    = (gamma ((v + d) / 2) * Dt) / ...
        ((pi * v)^(d / 2) * gamma(v / 2) * d_mahal);
end