close all
clearvars
clc


%% Caricare i dati e osservare la loro natura

% La generazione di numeri casuali avviene con algoritmi di diversa natura,
% qui impostiamo il generatore di default (Mersenne Twister) e poi
% impostiamo un seed pari a 1. Questo garantisce una omogeneita' nei
% risultati.

% Per casa e per curiosita' provate a ripetere due volte queste righe:
% rng('default')
% rng(1)
% randi(10,3)
% vedrete che la matrice 3x3 che avete generato sara' uguale.

rng('default');
rng(1);

load("data_lab8.mat")
figure(1)
plot(time,ecg)
xlabel('Time [s]')
ylabel('ECG [a.u.]')

%%  Calcolare la silhouette per ciasun K da 2 a 10 utilizzando le funzioni kmeans e silhouette (vedere l'help)
% Utilizzare il k-means clustering per suddividere i dati in gruppi (cluster)
% Impostare i parametri replicates a 10 (numero di volte in cui si deve ripetere 
% il clustering usando una nuova posizione iniziale del centroide)
% e il parametro maxiter con un massimo di 200 iterazioni. 
% Impostare come distanza la distanza euclidea.
% 
% Attenzione: la matrice in ingresso a k-means ha sulle righe i punti
% (soggetti) e sulle colonne le variabili


Klist=2:10;
S_mean=zeros(length(Klist),1);

for ii=1:length(Klist)

    % Seleziono il K i-esimo per il quale calcolare la silhouette
    K=Klist(ii);

    % k-means clustering
    [CLUST{ii}, Centroids{ii}] = kmeans(ecg, K,'replicates', 10,'maxiter',200,'distance','sqeuclidean');
    % CLUST: indici del clustering (vettore: size(ecg,1))
    % Centroids: locazioni dei centroidi dei vari clusters (matrice: k x size(ecg,2))
    
    figure(2)
    plot(time,Centroids{ii})
    drawnow
    title(['Centroidi per k=' num2str(K)])
    
    % Calcolo della silhouette 
    S = silhouette(ecg, CLUST{ii}); 

    S_mean(ii)= mean(S);

end

[max_S_mean,idx_max_S_mean]= max(S_mean);         
Kopt=Klist(idx_max_S_mean);

% Plottare il valore di silhouette media per ciascun K

figure(3)
plot(Klist, S_mean, '*-', 'LineWidth', 1 )
xlabel('Numero di clusters')
ylabel('Silhouette media')
hold on, plot(Kopt, max_S_mean, 'ro', 'LineWidth', 2) 
legend ('Silhouette media per differenti K', 'Silhouette media per Kopt')
title('Silhouette media al variare di K')

%% 3) Salviamo assegnazioni e centroidi del risutalto kmeans per k=kopt
% N.B. non ho risultati per k=1
CLUST_kmeans     = CLUST{Kopt-1};
Centroids_kmeans = Centroids{Kopt-1};

figure(4)
silhouette(ecg,CLUST_kmeans)
title(['K-means clustering con Kopt=' num2str(Kopt)])

% Rappresentare la matrice di similarità per il Kopt selezionato
% Ordinare i soggetti che appartengono ad ogni cluster (numero di clusters dato dal Kopt)
[~,id_ecg_sorted]=sort(CLUST_kmeans);
ecg_sorted = ecg(id_ecg_sorted,:);

% Calcolare la distanza euclidea tra ogni paio di osservazioni in ecg_sorted (p-dist)
Dist_vector= pdist(ecg_sorted);
% le distanze calcolate in Dist_vector sono nell'ordine (2,1), (3,1), ..., (m,1), (3,2), ..., (m,2), ..., (m,m–1).
% Dist_vector ha lunghezza m(m–1)/2, dove m è il numero di osservazioni

% Convertire Dist_vector in Dist_matrix (matrice mxm simmetrica con la
% diagonale principale a 0), la matrice delle distanze (squareform)
Dist_matrix=squareform(Dist_vector); 

% Mostrare in una figura la matrice di similarità (imagesc)
figure(5)
imagesc(1./Dist_matrix)
colorbar
colormap jet
title('Matrice di similarità')


%% Applicazione del Clustering Gerarchico

%%Dendrogramma dell'albero gerarchico
tree_ward = linkage(ecg,'ward');

figure(6)
dendrogram(tree_ward)
title('Hierarchical clustering tree (linkage: ward, no cutoff)')


figure(7)
cutoff = median([tree_ward(end-(Kopt-1),3) tree_ward(end-(Kopt-2),3)]);
dendrogram(tree_ward,'ColorThreshold',cutoff)
title('Hierarchical clustering tree (linkage: ward)')

tree_complete = linkage(ecg,'complete');
figure(8)
cutoff = median([tree_complete(end-(Kopt-1),3) tree_complete(end-(Kopt-2),3)]);
dendrogram(tree_complete,'ColorThreshold',cutoff)
title('Hierarchical clustering tree (linkage: complete)')

CC_complete = cophenet(tree_complete,pdist(ecg));
CC_ward     = cophenet(tree_ward,pdist(ecg));

%%

%%Ricavo l'assegnazione ai clusters con K=Kopt precedentemente selezionato
CLUST_tree = cluster(tree_ward,'maxclust',Kopt);
clear Centroids_tree

Centroids_tree = zeros(Kopt,size(ecg,2));
%%Calcolo i centroidi dei clusters identificati 
for k=1:Kopt
    idx_subj_clusterk= find(CLUST_tree==k);
    if length(idx_subj_clusterk)==1
        Centroids_tree(k,:)= ecg(idx_subj_clusterk,:);    
    else
        Centroids_tree(k,:)= mean(ecg(idx_subj_clusterk,:),'omitnan');
    end
end


%% Confronto dei risultati ottenuti dall%applicazione dei due metodi:
%%Calcolo la correlazione tra centroidi ottenuti con i due algoritmi di clustering
[CORR_centroids,p_centroids] = corr(Centroids_kmeans', Centroids_tree');

figure(9)
subplot(121)
imagesc(CORR_centroids)
colorbar
ylabel('Clustering K-means')
xlabel('Clustering Gerarchico')
title('Correlazione tra i centroidi')

subplot(122)
imagesc(p_centroids<0.05)
colorbar
colormap jet
ylabel('Clustering K-means')
xlabel('Clustering Gerarchico')
title('Significativita''')

%Calcolo le distanze tra i centroidi ottenuti con i due algoritmi di clustering
DIST_centroids = pdist2(Centroids_kmeans, Centroids_tree, 'euclidean');
figure(10)
imagesc(DIST_centroids)
colorbar
set(gca,'XTick',1:Kopt,'YTick',1:Kopt)
ylabel('Clustering K-means')
xlabel('Clustering Gerarchico')
title('Distanze tra centroidi')

[~,idx_centroids] = pdist2(Centroids_kmeans, Centroids_tree, 'euclidean','Smallest',1);

% Confrontare i centroidi e individuare quale centroide ricavato con il
% k-means si puo' ritenere simile o addirittura equivalente a quello
% ricavato con il clustering gerarchico

figure(11)
for c=1:Kopt
    subplot(3,2,c)
    plot(time,Centroids_kmeans(idx_centroids(c),:))
    hold on
    plot(time,Centroids_tree(c,:))
    title(['Gerarchico c=' num2str(c) ' Kmeans c=' num2str(idx_centroids(c))])
    xlabel('Time [s]')
    legend('K-means','Gerarchico')
end