function [mask P M S V] = tmm_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
%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 Kx1 matrix
global MU_INIT;

EPSILON = 1E-06;
ITER    = 20; %EM ITTERATIONS UPPER BOUND

N = size (x,1);

p1 = zeros (N,K); p_post = zeros (N,K); u = zeros (N,K);

P0 = zeros (K,1);
M0 = zeros (K,1);
S0 = zeros (K,1);
V0 = 5 + rand(K,1);

if isempty(MU_INIT)
    %INITIALIZE WITH K-MEANS
    mask = kmeans (x, K, 'emptyaction','singleton','display','off',...
        'start','cluster','replicates',20);
    for i=1:K
        ind = find (mask == i);
        if isempty (ind)
            P0(i) = eps;
            M0(i) = 0;
            S0(i) = 1000;
        else
            P0(i) = length(ind) / N;
            M0(i) = mean (x(ind));
            S0(i) = std (x(ind));
        end
    end
else
    M0 = MU_INIT;
    S0 = ones(K,1);
    P0 = (1/K) * ones (K,1);
end

%EM Loop
Qnew = inf;
for iter=1:ITER
    %Expectation
    for j=1:K
        p1(:,j) = t_pdf(x,M0(j),S0(j),V0(j))* P0(j);
    end
    px = sum (p1,2);
    px(find(px == 0)) = eps;
    for j=1:K
        p_post (:,j) = p1(:,j) ./ px;
        d_mahal = ((x - M0(j)).^2) / S0(j)^2;
        u (:,j) = (V0(j) + 1) ./ (V0(j) + d_mahal);
    end

    P = P0;
    M = M0;
    S = S0;
    V = V0;
    
    %convergence condition
    Qold = Qnew; Qnew = -sum(log(px));
    if Qold - Qnew < EPSILON
        break;
    end
    
    %Maximization
    P0  = sum(p_post,1)'/N;
    P0(find(P0 == 0)) = eps;
    for j=1:K
        M0(j) = sum(p_post(:,j).*u(:,j).*x) ...
            / sum (p_post(:,j).*u(:,j));
        S0(j) = sqrt (sum (p_post(:,j).*u(:,j).*...
            ((x- M0(j)).^2))/ sum (p_post(:,j)));
        if S0(j) < 1.0
            P0 (j) = eps;
            S0 (j) = 1000;
            continue;
        end
        V0(j) = get_dof (V0(j),1,u(:,j),p_post(:,j));
    end
end


% normalization of priors
P = P / sum(P);

for j=1:K
    p1(:,j) = t_pdf(x,M(j),S(j),V(j))* P(j);
end
px = sum (p1,2);
px (find(px == 0)) = eps;
L = sum (px);

for j=1:K
    p_post (:,j) = p1(:,j) ./ px;
end
[q mask] = max (p_post,[],2);

% -----

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

EPSILON = .001; MAX_ITER = 50;

optim_par = optimset ('TolX',1E-04,'TolFun',1E-06,...
    'Algorithm','active-set',...
    'LargeScale','Off','GradObj','Off',...
    'MaxFunEvals',5000,...
    'Display','Off');

u(find(u<=0)) = eps;
n = sum(t); Q = -psi(v/2)+log(v/2)+1+sum(t.*(log(u)-u))/n;
fun_z = @(x)((psi((x+p)/2) - log((x+p)/2) + Q)^2);

%ROOT FINDING => GLOBAL MINIMIZATION OF FUNCTION^2
%ALTERNATIVE WE CAN USE FZERO BUT!!! WE MUST GUARANTEE
%THAT POSITIVE Xs ARE SEARCHED

fval = inf;
for iter=1:MAX_ITER
    x0 = 1 + 8*rand(1,1);
    [new_v1 fval1] = fmincon (fun_z,x0,[],[],[],[],EPSILON,10,[],optim_par);
    
    if fval1 < fval
        fval = fval1;
        new_v = new_v1;
    end
    if fval < EPSILON break; end
end

%
% function: t_pdf
%
function [tp] = t_pdf (x, m, s, v)
d_mahal = ((x-m)/s).^2;
tp = gamma((v+1)/2) ./ (s*sqrt(pi*v)*gamma(v/2)*...
    (1+d_mahal/v).^((v+1)/2));