%% LABORATORIO 11 - METODI STATISTICI PER LA BIOINGEGNERIA
%% Canale B - A.A. 2024/2025
%% Docente: Enrico Longato

clc
clear all
close all

% Inizializzazione RNG
rng default
rng(42)

%% Caricamento dati
load data_lab11


%% Visualizzazione delle tracce ECG
signals_to_plot = [1 14 23];

figure
subplot(211)
plot(time, ecg) % MATLAB intuisce cosa fare dalla dimensione delle matrici
xlabel('Tempo (s)')
ylabel('ECG (a. u.)') % Arbitrary units (nel senso che non le sappiamo)
subtitle('Tutti i segnali')

subplot(212)
hold on
for i_signal = signals_to_plot % Ricordo che il ciclo for può scandire vettori arbitrari
    plot(time, ecg(i_signal, :), 'DisplayName', num2str(i_signal)) % DisplayName serve per la legenda
end
xlabel('Tempo (s)')
ylabel('ECG (a. u.)') 
subtitle('Segnali selezionati')
legend % Attiva la legenda

%% K-means con 3 cluster e silhouette
K_fixed = 3; % Fisso K = 3 numero di cluster
N_replicates = 10; % Ripetizioni del K-means per non finire in ottimi locali
N_iter = 200; % Numero massimo di iterazioni dell'algoritmo per ogni replica

[cluster_fixed, centroid_fixed] = kmeans(ecg, K_fixed, ...
                                         'Distance', 'sqeuclidean', ...
                                         'Replicates', N_replicates, ...
                                         'MaxIter', N_iter);

% Plot dei centroidi
figure
hold on
for i_centroid = 1:K_fixed % Ricordo che il ciclo for può scandire vettori arbitrari
    plot(time, centroid_fixed(i_centroid, :), ...
        'DisplayName', num2str(i_centroid)) % DisplayName serve per la legenda
    ylim([-1 2])
end
xlabel('Tempo (s)')
ylabel('ECG (a. u.)')
title(['Centroidi con K = ' num2str(K_fixed)])
legend % Attiva la legenda

% Calcolo della silhouette media
sil_fixed = silhouette(ecg, cluster_fixed); % Silhouette di tutti i punti (cioè di tutte le tracce ECG)
sil_mean_fixed = mean(sil_fixed); % Faccio la media

% Si può anche rappresentare graficamente togliendo gli argomenti di uscita
figure
silhouette(ecg, cluster_fixed)
title(['Silhouette di ciascun segnale (media = ', num2str(sil_mean_fixed), ')'])

%% K-means alla ricerca del numero di cluster ottimo tra 2 e 10 inclusi
% Si tratta di fare un ciclo for copiando le istruzioni di prima e
% salvandole opportunamente a ogni iterazione e, poi, fare il plot dei
% centroidi corrispondenti al cluster ottimo (che vanno salvati in un cell
% array)

K_possible = 2:10; % Numero di cluster da scandire
silhouette_vector = []; % Vettore vuoto delle silhouette medie

% Cell array dei centroidi
centroid_cell_array = {}; % Dentro a ciascuna cella ci sarà una matrice [k_curr x 300]

for k_curr = K_possible

    % Applicare l'algoritmo K-means con i parametri di prima, ma col numero
    % di cluster corrente
    [cluster_curr, centroid_curr] = kmeans(ecg, k_curr, ...
                                           'Distance', 'sqeuclidean', ...
                                           'Replicates', N_replicates, ...
                                           'MaxIter', N_iter);
    
    % Calcolo i valori di silhouette corrispondenti
    sil_curr = silhouette(ecg, cluster_curr);
    sil_mean_curr = mean(sil_curr);

    % Memorizzo il necessario
    silhouette_vector = [silhouette_vector sil_mean_curr]; % Accrescimento dinamico della lista
    centroid_cell_array = [centroid_cell_array centroid_curr]; % Accrescimento dinamico del cell array

    % NB: potevo memorizzare anche sil_curr in un cell array, ma non era
    % richiesto il plot finale

end

% Trovo il massimo della silhouette e rappresento graficamente
[max_sil, i_max_sil] = max(silhouette_vector); % Attenzione! i è l'indice che parte da 1
k_max_sil = K_possible(i_max_sil); % Attenzione! Per il plot devo trovare il valore di K corrispondente

figure
hold on
plot(K_possible, silhouette_vector, 'k', 'Marker', '.')
plot(k_max_sil, max_sil, 'ro')
xlabel('Numero di cluster')
ylabel('Silhouette media')
title('Silhouette media al variare del numero di cluster')

% Disegno i centroidi
figure
plot(time, centroid_cell_array{i_max_sil}) % Sempre all'indice massimo; è l'altro modo (con la matrice) di fare il plot
xlabel('Tempo (s)')
ylabel('ECG (a. u.)')
title(['Centroidi con K = ' num2str(k_max_sil)])
ylim([-1 2])