% Reinizializzare a default il generatore dei numeri casuali
% Generare una matrice di dimension n1-by-p (n1=10 e p=2) di numeri casuali
% dalla distribuzione N(0,1). Generare una matrice di dimension n2-times-p
% (n2=20 e p=2) di numeri casuali dalla distribuzione N(0+dx,1) con dx=5.
%
% Impostare a 3 il numero dei gruppi. Utilizzare come valore iniziale dei
% centroidi le unità 5, 15 e 25. Calcolare le distanze di ogni unita da
% questi 3 centroidi. Assegnare ogni unità al gruppo con la distanza
% Euclidea più piccola. Trovare i nuovi valori dei centroidi e calcolare il
% valore della funzione obiettivo del metodo delle k-medie (minimizzazione
% della devianza all'interno dei gruppi).
%
% Effettuare 10 iterazioni e monitorare il valore della funzione obiettivo.
% Rappresentare graficamente il valore della funzione obiettivo al variare
% delle iterazioni
%
% Applicare la funzione kmeans con k=3 e numero di repliche pari a 100
% Mostrare la classificazione finale ed i centroidi


rng default
% Generazione dei dati
n1=10;
n2=20;
p=2;
X1=randn(n1,p);
dx=5;
X2=dx+randn(n2,p);
% Combinare i due gruppi in una matrice di dimensione (n1+n2)-by-p;
X=[X1; X2];

%% Una sola iterazione
% index = vettore che contiene gli indici da utilizzare come valori iniziali dei centroidi
index = [5 15 25];

k=3;
n=size(X,1);

% cini= matrice di dimensione 3xp che contiene
% i centroidi iniziali dei 3 gruppi
cini=X(index,:);

% D = matrice di dimensione nxk che conterrà le distanze di ciascuna unità
% da ciascuno dei k centroidi
% Le righe di D si riferiscono alle unità
% le colonne di D si riferiscono ai clusters
D=zeros(n,k);


% Viene calcolata la distanza (al quadrato) di ogni unità dai k centroidi
% D(i,j) contiene la distanza Euclidea dell'unità i dal cluster j
% D(i,:) (vettore 1 x k) che contiene le distanze dell'unità i dai k centroidi
% D(:,j) (vettore n x 1) che contiene le distanze delle n unità dal cluster j
for j=1:k
    D(:,j)=sum((X-cini(j,:)).^2,2);
end

% dist = vettore di dimensione n x 1 che contiene le distanze di ogni unità
% dal centroide più vicino. ind = vettore di dimensione n x 1 che contiene
% gli assegnamenti
[dist,ind]=min(D,[],2);

% Ricalcolo dei nuovi k centroidi e determino il valore della funzione
% obiettivo
% sumd = valore della funzione obiettivo (da minimizzare)
% Somma delle distanze (al quadrato) di ciascuna unità dal rispettivo centroide
sumd=0;
for j=1:k
    Xj=X(ind==j,:);
    cini(j,:)=mean(Xj);
    sumd=sumd+sum(sum((Xj-cini(j,:)).^2,2));
end

disp('Valore della funzione obiettivo alla prima iterazione')
% sum è il valore della devianza nei gruppi
disp(sumd)
%% 10 iterazioni
% index = vettore che contiene gli indici da utilizzare come valori iniziali dei centroidi
index = [5 15 25];

k=length(index);
n=size(X,1);

cini=X(index,:);

D=zeros(n,k);

niter=10;
% Obj = vettore che contiene il valore della funzione obiettivo per ogni
% iterazione
Obj=zeros(niter,1);

for iter =1:niter

    for j=1:k
        D(:,j)=sum((X-cini(j,:)).^2,2);
    end

    [dist,ind]=min(D,[],2);

    % Ricalcolo dei nuovi centroidi e del nuovo valore della funzione obiettivo
    sumd=0;
    for j=1:k
        Xj=X(ind==j,:);
        cini(j,:)=mean(Xj);
        sumd=sumd+sum(sum((Xj-cini(j,:)).^2,2));
    end

    % Il valore della funzione obiettivo all'iterazione iter viene inserita
    % in posizione iter del vettore Obj
    Obj(iter)=sumd;

end

plot(Obj)
xlabel('Iterazioni')
ylabel('Valore della funzione obiettivo')
% print -depsc figs\kmeansINI.eps;


%% Utilizzo funzione kmeans
k=3;
[idx,cent,sumdj]=kmeans(X,k,"Replicates",100);
% Verifico che il secondo argomento di output della funzione kmeans
% contiene le medie aritmetiche dei 3 gruppi
Centroidi=grpstats(X,idx);
assert(max(abs(cent-Centroidi),[],"all")<1e-08,'Errore di programmazione')
sumd=sum(sumdj);
disp(['Il valore della funzione obiettivo minimizzata è=' num2str(sumd)])

disp("Centroidi dei gruppi")
disp(cent)

% Classificazione finale e centroidi
gscatter(X(:,1),X(:,2),idx,'brk','xos',15)
hold('on')
gscatter(Centroidi(:,1),Centroidi(:,2),1:k,'brk','p',15)
title('Classificazione finale e centroidi')
% print -depsc figs\kmeansINI1.eps;