% Lezione di lab TTR  -27/11/2013- 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)];
 
figure;
scatter(v(:,1),v(:,2)); hold on

%[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;

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...





