% Lezione di lab TTR  -18/04/2016- Expectation Maximization - Marco Cristani


clear all;

K=2;

v = [mvnrnd([0,0]',[1, 0; 0,1],500);...
     mvnrnd([2,3]',[2, 0; 0,1],500)];
 


%[label, model, llh] = emgm(v', K);
[model, llh] = emgm(v', K);
figure;
plot(llh); set(gcf,'name','EM - log likelihood evolution')
xlabel('Iterations')
ylabel('loglikelihood')

level = 2.5;
dashed = 0;

figure;
scatter(v(:,1),v(:,2)); hold on
for k=1:K
y = elipsnorm(model.mu(:,k),model.Sigma(:,:,k),level,dashed);
end

%% ESERCIZIO 1: creare due classi di misture di gaussiane, addestrarle, fare prove di classificazione

clear all; 

% Prima classe;
pi1_tmp=0.5;
pi1=[pi1_tmp,1-pi1_tmp];
mu1= [[0,0]',[1,1]'];
Sigma1(:,:,1)= [1,0;0,1];
Sigma1(:,:,2)= [0.1,0;0,0.5];

% Seconda classe;
pi2_tmp = 0.2;
pi2=[pi2_tmp,1-pi2_tmp];
mu2= [[2,2]',[3,1]'];
Sigma2(:,:,1)= [2,0;0,1];
Sigma2(:,:,2)= [0.5,0;0,0.1];


% Generazione campionidi train
n_p= 500; % numero punti per ogni classe
data1=[];
num_c1=0;
num_c2=0;
for i=1:n_p
    seed = unifrnd(0,1);
    if seed > pi1_tmp
        data1 = [data1, mvnrnd(mu1(:,1),Sigma1(:,:,1))'];
        num_c1=num_c1+1;
    else
        data1 = [data1, mvnrnd(mu1(:,2),Sigma1(:,:,2))'];
        num_c2=num_c2+1;
    end
end

data2=[];
for i=1:n_p
    seed = unifrnd(0,1);
    if seed > pi2_tmp
        data2 = [data2, mvnrnd(mu2(:,1),Sigma2(:,:,1))'];
    else
        data2 = [data2, mvnrnd(mu2(:,2),Sigma2(:,:,2))'];
    end
end

data = [data1,data2];
labels = [ones(1,n_p),ones(1,n_p).*2];

figure

set(gcf,'name','Dati iniziali');
scatter(data(1,:),data(2,:),[],labels);
hold on

% ---Addestramento---
K=2;
level = 2.5;


%classe 1
[model1, llh1] = emgm(data1, K);
dashed = 0;
for k=1:K
y = elipsnorm(model1.mu(:,k),model1.Sigma(:,:,k),level,dashed);
end

%classe 2
[model2, llh2] = emgm(data2, K);
dashed = 1;
for k=1:K
y = elipsnorm(model2.mu(:,k),model2.Sigma(:,:,k),level,dashed);
end


% Generazione campionidi test
n_p= 50; % numero punti per ogni classe
data1t=[];
num_c1=0;
num_c2=0;
for i=1:n_p
    seed = unifrnd(0,1);
    if seed > pi1_tmp
        data1t = [data1t, mvnrnd(mu1(:,1),Sigma1(:,:,1))'];
        num_c1=num_c1+1;
    else
        data1t = [data1t, mvnrnd(mu1(:,2),Sigma1(:,:,2))'];
        num_c2=num_c2+1;
    end
end

data2t=[];
for i=1:n_p
    seed = unifrnd(0,1);
    if seed > pi2_tmp
        data2t = [data2t, mvnrnd(mu2(:,1),Sigma2(:,:,1))'];
    else
        data2t = [data2t, mvnrnd(mu2(:,2),Sigma2(:,:,2))'];
    end
end

datat = [data1t,data2t];
labelst = [ones(1,n_p),ones(1,n_p).*2];

set(gcf,'name','Dati iniziali');
scatter(datat(1,:),datat(2,:),[],labelst,'Marker','x');


% --- Testing ---

for i=1:size(datat,2)
    for k=1:K
        l1tmp(k) = model1.weight(k).*mvnpdf(datat(:,i),...
            model1.mu(:,k),model1.Sigma(:,:,k));
    end
    ll1(i)= log(sum(l1tmp));
    
    for k=1:K
        l2tmp(k) = model2.weight(k).*mvnpdf(datat(:,i),...
            model2.mu(:,k),model2.Sigma(:,:,k));
    end
    ll2(i)= log(sum(l2tmp));
end

[thrash,labelsMog]= max([ll1;ll2]);

% --- Valutazione ---
conf=zeros(2,2);

for i=1:size(datat,2)
conf(labelsMog(i),labelst(i))= conf(labelsMog(i),labelst(i))+1;
end



%% ESERCIZIO 2: applicare EM alle facce, come metodologia di clustering

clear all;
dire    =    '.\images'; % ORL face database
list    =   dir(strcat(dire,'\*.bmp'));
M       =   size(list,1);
tmp     =   imread(strcat(dire,'\',list(1).name));
[r,c]   =   size(tmp);
TMP     =   zeros(r*c,M);
figure;
for i=1:M
    tmp = imread(strcat(dire,'\',list(i).name));
%     imagesc(tmp); axis image; colormap gray;
%     drawnow;
%     pause(0.01);
    TMP(:,i)= tmp(:);
%     mont(:,:,1,i)= tmp;
end
% 
% figure;
% montage(mont);
% set(gcf,'name','Gallery');
TMP                     =   double(TMP);
media                   =   mean(TMP,2);
A(:,:)                  =   TMP-repmat(media,1,M);

[U,lambda]              =   eigen_training(A);


num_eig           =   40;
omega       =   U(:,1:num_eig)'*A; % proiezione;
figure; set(gcf,'name','Visualizzazione 2D dei volti');
scatter(omega(1,:),omega(2,:));

l           =   reshape(repmat([1:40],10,1),400,1);
for i=1:M
    text(omega(1,i),omega(2,i),num2str(l(i)))
end
hold on

K = 5;
[model, llh] = emgm(omega, K);

level = 2.5;
dashed = 0;

for k=1:K
y = elipsnorm(model.mu(1:2,k),model.Sigma(1:2,1:2,k),level,dashed);
end

%% ESERCIZIO 3: visualizzare le medie

for k=1:K
rec(:,k) = media + U(:,1:num_eig)*model.mu(:,k);
end

figure; set(gcf,'name','Medie dei Cluster');
num_plot= ceil(sqrt(K));
for k=1:K
subplot(num_plot,num_plot,k);
imagesc(reshape(rec(:,k),r,c)); axis image; colormap gray;
end

%% ESERCIZIO 4: visualizzare i campioni appartenenti a ciascun cluster
%bisogna modificare il codice dell'EM...


%% ESERCIZIO 5: eseguire prove di classificazione su biondi, mori e calvi
clear all;

labelm_tmp = [1,21,31,51,61,71,81,91,101, 131,161,181,191,201,211,221,231,241,251,281,291,301,321,351,371,381];
labelb_tmp=[11,41,111,141,171,311,341,391];

labelc_tmp=[121,151,261,271,331,361];

label =[]
for i=1:length(labelm_tmp)
    label(labelm_tmp(i):labelm_tmp(i)+9)=1;
end
for i=1:length(labelb_tmp)
    label(labelb_tmp(i):labelb_tmp(i)+9)=2;
end
for i=1:length(labelc_tmp)
    label(labelc_tmp(i):labelc_tmp(i)+9)=3;
end

dire    =    '.\images'; % ORL face database
list    =   dir(strcat(dire,'\*.bmp'));
M       =   size(list,1);
tmp     =   imread(strcat(dire,'\',list(1).name));
[r,c]   =   size(tmp);
TMP     =   zeros(r*c,M);
for i=1:M
    tmp = imread(strcat(dire,'\',list(i).name));
    TMP(:,i)= tmp(:);
end
TMP                     =   double(TMP);
media                   =   mean(TMP,2);
A(:,:)                  =   TMP-repmat(media,1,M);

[U,lambda]              =   eigen_training(A);

T           =   300;
X           =   U(:,1:T)'*A; % proiezione;

[d,N] = size(X);
K = max(label); % numero classi in gioco;

% 1. determino le classi Ck
for k = 1:K
    a = find (label == k);
    Ck{k} = X(:,a);
end
% 2. determino le medie
for k = 1:K
    mk{k} = mean(Ck{k},2);
end
% 3. determino la numerosità della classe
for k = 1:K
    [d, Nk(k)] = size(Ck{k});
end
% 4. determino le within class scatter
for k = 1:K
    S{k} = 0;
    for i = 1:Nk(k)
        S{k} = S{k} + (Ck{k}(:,i)-mk{k})*(Ck{k}(:,i)-mk{k})';
    end
    S{k} = S{k}./Nk(k);
end
Swx = 0;
for k = 1:K
    Swx = Swx + S{k};
end

% 5. determino la between class covariance
% 5.1 determino la media totale
m = mean(X,2);
Sbx = 0;
for k=1:K
    Sbx = Sbx + Nk(k)*((mk{k} - m)*(mk{k} - m)');
end
Sbx = Sbx/K;

MA = inv(Swx)*Sbx;

% eigenvalues/eigenvectors
[V,D] = eig(MA);

% 5: transform matrix
A = V(:,1:3);

% 6: transformation
Y = A'*X;

% 7: plot
figure, scatter(Y(1,:),Y(2,:),[],label)
for i=1:M
    text(Y(1,i),Y(2,i),num2str(label(i)))
end


