load cities;
nRow = 329;
inputs = [ratings(:,1), ratings(:,2);];

C = cov(ratings);
C = corr(ratings,ratings);

f1 = figure;
clima = ratings(1:329,1);
figClima = histogram(clima);
legend('Clima');

f2 = figure;

abitazioni = ratings(1:329,2);
figabitazioni = histogram(abitazioni);
legend('abitazioni');

f3 = figure;
salute = ratings(1:329,3);
figsalute = histogram(salute);
legend('salute');

f4 = figure;
criminalita = ratings(1:329,4);
fig_criminalita = histogram(criminalita);
legend('criminalita');

f5 = figure;
trasporti = ratings(1:329,5);
fig_trasporti = histogram(trasporti);
legend('trasporti');

f6 = figure;
istruzione = ratings(1:329,6);
fig_istruzione = histogram(istruzione);
legend('istruzione');

f7 = figure;
arte = ratings(1:329,7);
fig_arte = histogram(arte);
legend('arte');

f8 = figure;
svago = ratings(1:329,8);
fig_svago = histogram(svago);
legend('svago');

f9 = figure;
economia = ratings(1:329,9);
fig_economia = histogram(economia);
legend('economia');


Setto il numero di cluster
k = 4;



disp('--------------------');
[climaMin,indexMinClima] = min(clima);
[climaMax,indexMaxClima] = max(clima);

cityWorstClima = names(indexMinClima,:);
cityBestClima = names(indexMaxClima,:);

disp(['city Worst Clima:', cityWorstClima]);
disp(['city Best Clima:', cityBestClima]);

[abitazioniMin,indexMinabitazioni] = min(abitazioni);
[abitazioniMax,indexMaxabitazioni] = max(abitazioni);

disp('--------------------');
cityWorstabitazioni = names(indexMinabitazioni,:);
cityBestabitazioni = names(indexMaxabitazioni,:);

disp(['city Worst abitazioni:', cityWorstabitazioni]);
disp(['city Best abitazioni:', cityBestabitazioni]);

disp('--------------------');

[saluteMin,indexMinsalute] = min(salute);
[saluteMax,indexMaxsalute] = max(salute);

cityWorstsalute = names(indexMinsalute,:);
cityBestsalute = names(indexMaxsalute,:);

disp(['city Worst salute:', cityWorstsalute]);
disp(['city Best salute:', cityBestsalute]);

[criminalitaMin,indexMincriminalita] = min(criminalita);
[criminalitaMax,indexMaxcriminalita] = max(criminalita);

cityWorstcriminalita = names(indexMincriminalita,:);
cityBestcriminalita = names(indexMaxcriminalita,:);

disp('--------------------');
disp(['city Worst criminalita:', cityWorstcriminalita]);
disp(['city Best criminalita:', cityBestcriminalita]);

[trasportiMin,indexMintrasporti] = min(trasporti);
[trasportiMax,indexMaxtrasporti] = max(trasporti);

cityWorsttrasporti = names(indexMintrasporti,:);
cityBesttrasporti = names(indexMaxtrasporti,:);

disp('--------------------');
disp(['city Worst trasporti:', cityWorsttrasporti]);
disp(['city Best trasporti:', cityBesttrasporti]);


[istruzioneMin,indexMinistruzione] = min(istruzione);
[istruzioneMax,indexMaxistruzione] = max(istruzione);

cityWorstistruzione = names(indexMinistruzione,:);
cityBestistruzione = names(indexMaxistruzione,:);

disp('--------------------');
disp(['city Worst istruzione:', cityWorstistruzione]);
disp(['city Best istruzione:', cityBestistruzione]);

[arteMin,indexMinarte] = min(arte);
[arteMax,indexMaxarte] = max(arte);

cityWorstarte = names(indexMinarte,:);
cityBestarte = names(indexMaxarte,:);

disp('--------------------');
disp(['city Worst arte:', cityWorstarte]);
disp(['city Best arte:', cityBestarte]);

[svagoMin,indexMinsvago] = min(svago);
[svagoMax,indexMaxsvago] = max(svago);

cityWorstsvago = names(indexMinsvago,:);
cityBestsvago = names(indexMaxsvago,:);

disp('--------------------');
disp(['city Worst svago:', cityWorstsvago]);
disp(['city Best svago:', cityBestsvago]);

[economiaMin,indexMineconomia] = min(economia);
[economiaMax,indexMaxeconomia] = max(economia);

cityWorsteconomia = names(indexMineconomia,:);
cityBesteconomia = names(indexMaxeconomia,:);

disp('--------------------');
disp(['city Worst economia:', cityWorsteconomia]);
disp(['city Best economia:', cityBesteconomia]);

t = [ratings(1:20,1) ratings(1:20,2)];

scatter(t(:,1) , t(:,2));

M = mean(ratings);

S = std(ratings);

D = cov(ratings);

% ---------- CALCOLO MANUALE VARIANZA  ------------
% citta da prendere, per l'esempio solo 5
valueToTake = 5;
% Prendo la prima colonna, fino al 5 valore
% Prendo la seconda colonna fino al 5 valore
% Infine combino i vettori in un unica matrice 5x2
A = [ratings(1:valueToTake,1) ratings(1:valueToTake,2)];

% Calcolo la media (per colonne)
M = mean(A);
% Calcolo la varianza (per colonna)
% S = std(A); Non serve
% Creo un vettore di tutti zeri (con 5 righe)
X1 = zeros(size(ratings(1:valueToTake,1)));
% Combino questo vettore vuoto con la matrice per avere uno spazio per i
% calcoli
B = [A X1 X1 X1];

% Calcolo (valore 'ennesimo' - media) e lo salvo nella sua colonna
for i = 1:valueToTake
        B(i,3) = B(i,1)-M(1,1);
end

t = M(1,2);

% Calcolo (valore 'ennesimo' - media) e lo salvo nella sua colonna
for i = 1:valueToTake
        B(i,4) = B(i,2)-M(1,2);
end

% Calcolo la moltiplicazione tra i 2 fattori e lo salvo nella sua colonna
for i = 1:valueToTake
        B(i,5) = (B(i,3)) *(B(i,4));
end

% Calcolo covarianza
totColonna = 0;
for i = 1:valueToTake
        t = (B(i,5));
        totColonna = totColonna +t;
end

CovarianzaAB = totColonna / (valueToTake-1);

% Matlab
x1 = ratings(1:valueToTake,1);
x2 = ratings(1:valueToTake,2);
D = cov(x1,x2);

% ---------- FINE CALCOLO MANUALE VARIANZA  ------------



% % ---------- CALCOLO MANUALE MATRICE CORRELAZIONE  ------------

% citta da prendere, per l'esempio solo 5
valueToTake = 329;

% Prendo la prima colonna, fino al 5 valore
% Prendo la seconda colonna fino al 5 valore
% Infine combino i vettori in un unica matrice 5x2
x1 = ratings(1:valueToTake,1);
x2 = ratings(1:valueToTake,2);
A = [x1 x2];
% stampo il grafico
% scatter(x1,x2);
% Calcolo la media delle 2 features
Mvector = mean(A);
M = [Mvector(1,1); Mvector(1,2)];

% Creo l'oggetto cell dove mi posso salvare i valori intermedi dei calcoli
CellArrey = cell(valueToTake, 3);

%Calcolo vero a proprio
for i = 1:valueToTake

    % 1) creo la matrice (2x1) con i 2 valori
    temp_mat  = [x1(i,1); x2(i,1)];
    CellArrey{i,1} = temp_mat;

    % 2) salvo la differenza tra il valore del punto attuale e la media
    CellArrey{i,2} = minus(temp_mat,M);
    G = minus(temp_mat,M);

    %3) Mi calcolo la matrice appena trovata per se stessa trasposta
    Gtrasposta = G.';
    Step3 = G * Gtrasposta;
    CellArrey{i,3} = Step3;

end

% Creo la matrice per il calcolo finale
CovarianzaMatrix = [0,0;0,0];

% Calcolo la matrice di covarianza, sommando tutti i valori parziali e
% dividendo il risultato per n-1

% Calcolo la matrice di covarianza, sommando tutti i valori parziali
for i = 1:valueToTake
    % Mi salvo temporaneamente la matrice salvata nell'oggetto CellArrey
    tempCell = CellArrey{i,3};
    % Mi salvo la somma parziale progressiva di tutti i valori:
    % Riga=1 e Colonna=1
    CovarianzaMatrix(1,1) = CovarianzaMatrix(1,1) + tempCell(1,1);
    % Riga=1 e Colonna=2
    CovarianzaMatrix(1,2) = CovarianzaMatrix(1,2) + tempCell(1,2);
    % Riga=1 e Colonna=2
    CovarianzaMatrix(2,1) = CovarianzaMatrix(2,1) + tempCell(2,1);
    % Riga=2 e Colonna=2
    CovarianzaMatrix(2,2) = CovarianzaMatrix(2,2) + tempCell(2,2);
end

% dividendo il risultato per n-1
CovarianzaMatrix(1,1) = CovarianzaMatrix(1,1) / (valueToTake-1);
CovarianzaMatrix(1,2) = CovarianzaMatrix(1,2) / (valueToTake-1);
CovarianzaMatrix(2,1) = CovarianzaMatrix(2,1) / (valueToTake-1);
CovarianzaMatrix(2,2) = CovarianzaMatrix(2,2) / (valueToTake-1);

MatriceCorrelazione = [0,0;0,0];

% Calcolo la matrice di correlazione
% Si prendono gli indici i e j e si calcola per ogni cella il valore
% con la formula (ij)/((ii)*(jj))^1/2
for i = 1:2
    for j = 1:2
        MatriceCorrelazione(i,j) = CovarianzaMatrix(i,j) / sqrt(CovarianzaMatrix(i,i) * CovarianzaMatrix(j,j));
    end
end

% % ---------- FINE CALCOLO MANUALE MATRICE CORRELAZIONE  ------------

% Matlab
E = corr(ratings,ratings);


% % ---------- K MEANS  ------------

% Divido i dati in k classi
KMeansResult = kmeans(ratings,k);

% Divido i dati di ogni features in k classi
% La label del cluster segue un ordine crescente
Kclima = orderCluster(clima,k,nRow);
Kabitazioni = orderCluster(abitazioni,k,nRow);
Ksalute = orderCluster(salute,k,nRow);
Kcriminalita = orderCluster(criminalita,k,nRow);
Ktrasporti = orderCluster(trasporti,k,nRow);
Kistruzione = orderCluster(istruzione,k,nRow);
Karte = orderCluster(arte,k,nRow);
Ksvago = orderCluster(svago,k,nRow);
Keconomia = orderCluster(economia,k,nRow);

% Compongo una matrice con tutte le features
KRes = [Kclima Kabitazioni Ksalute Kcriminalita Ktrasporti Kistruzione Karte Ksvago Keconomia];

% mean - mode

% Scelgo come etichetta del cluster la colonna con frequenza maggiore
ModeKRes = mode(KRes,2);

% Calcolo una colonna con la media
MediaKRes = mean(KRes,2);
% Arrotondo alla cifra più vicina
for i = 1:329
        x = MediaKRes(i,1);
        MediaKRes(i,1) = round(x);
end

% Conto quante somiglianze ci sono con il Kmean
sameValue = 0;
compareVector = zeros(size(ratings(1:329,1)));

for i = 1:329
    valueKmeans = KMeansResult(i,1);
    valueMeanSingleKmeans = MediaKRes(i,1);
    if valueKmeans == valueMeanSingleKmeans
            sameValue = sameValue +1;
            compareVector(i,1) = 1;
    end
end
Kcompare = [KMeansResult MediaKRes compareVector];
% Conto quante volte i valori sono stati classificati nello stesso cluster
similarity = (sameValue/329)*100;

fig_KMeansResult = histogram(KMeansResult);
legend('KMeans');
f1 = figure;

fig_MediaKRes = histogram(MediaKRes);
legend('Media dei KMean');
f2 = figure;

% creates a Davies-Bouldin criterion clustering evaluation object.

eva1 = evalclusters(ratings,KMeansResult,'DaviesBouldin');
eva2 = evalclusters(ratings,MediaKRes,'DaviesBouldin');
eva3 = evalclusters(ratings,ModeKRes,'DaviesBouldin');

eva4 = evalclusters(ratings,'kmeans','DaviesBouldin','klist',[1:6]);

% ---------- FINE K MEANS  ------------

% % ---------- KMEANS CONFRONTO 2 FEATURES  ------------

% Decido il K
k = 2;
% Creo una matrice con solo 2 features
twoFeatures =[ratings(:,1), ratings(:,2)];
% Computo il kmeans
KMeansResult = kmeans(twoFeatures,k);
% Stampo il risultato
f5 = figure;1
gscatter(twoFeatures(:,1),twoFeatures(:,2),KMeansResult,'rbg','xod');

% Stima della bonta dei cluster con 'DaviesBouldin'
E = evalclusters(ratings,'kmeans','DaviesBouldin','klist',[1:6]);
% Stampa l'errore
f3 = figure;
plot(E);
% Stampa i dati colorandoli in base al cluster di appartenenza
f4 = figure;
gscatter(twoFeatures(:,1),twoFeatures(:,2),E.OptimalY,'rbg','xod');

% % ---------- FINE KMEANS CONFRONTO 2 FEATURES  ------------



% ---------- Fuzzy c-means CONFRONTO 2 FEATURES  ------------

% Decido il K
k = 2;
% Creo una matrice con solo 2 features
fcmdata =[ratings(:,1), ratings(:,2)];

t = size(fcmdata);
nRow = t(1,1);
% Mi creo un vettore di appoggio con tutti 0
finalVector = zeros(nRow,1);

% Computo il Fuzzy c-means
[centers,U] = fcm(fcmdata,k);

maxU = max(U);
index1 = find(U(1,:) == maxU);
index2 = find(U(2,:) == maxU);

plot(fcmdata(index1,1),fcmdata(index1,2),'ob')
hold on
plot(fcmdata(index2,1),fcmdata(index2,2),'or')
plot(centers(1,1),centers(1,2),'xb','MarkerSize',15,'LineWidth',3)
plot(centers(2,1),centers(2,2),'xr','MarkerSize',15,'LineWidth',3)
hold off

index1 = index1';
t = size(index1);
countLabel = t(1,1);
for i = 1:countLabel
    finalVector(index1(i,1),1)=1;
end

index2 = index2';
t = size(index2);
countLabel = t(1,1);
for i = 1:countLabel
    finalVector(index2(i,1),1)=2;
end
% Stima della bonta dei cluster con 'DaviesBouldin'
eva5 = evalclusters(fcmdata,finalVector,'DaviesBouldin');

% ---------- FINE Fuzzy c-means CONFRONTO 2 FEATURES  ------------

% ---------- CONFRONTO TRA I 2 ALGORITMI ------------

table = zeros(10,3);

for k = 2:10
    % Kmeans
    KMeansResult = kmeans(ratings,k);
    eva1 = evalclusters(ratings,KMeansResult,'DaviesBouldin');

    % Fuzzy c-means
    [centers,U] = fcm(ratings,k);
    maxU = max(U);
    finalVector = zeros(nRow,1);

    for i = 1:k
     index = find(U(i,:) == maxU);
     index = index';
     t = size(index);
     countLabel = t(1,1);

     for j = 1:countLabel
         finalVector(index(j,1),1)=i;
     end

    end

    eva2 = evalclusters(ratings,finalVector,'DaviesBouldin');

    table(k,1) = k;
    table(k,2) = eva1.CriterionValues;
    table(k,3) = eva2.CriterionValues;
end


% ---------- Self-Organizing Map ------------

% Create a Self-Organizing Map
dimension1 = 2;
dimension2 = 2;
net = selforgmap([dimension1 dimension2]);

% Train the Network
[net,tr] = train(net,inputs);

% Test the Network
outputs = net(inputs);

% View the Network
view(net)

% Plots
% Uncomment these lines to enable various plots.
% figure, plotsomtop(net)
% figure, plotsomnc(net)
% figure, plotsomnd(net)
% figure, plotsomplanes(net)
% figure, plotsomhits(net,inputs)
% figure, plotsompos(net,inputs)
%
%
% save('ratings.mat','ratings');
%
% Solve a Clustering Problem with a Self-Organizing Map
% Script generated by NCTOOL
%
% This script assumes these variables are defined:

%   simpleclusterInputs - input data.