% Klasifikace pomoci dat z velkeho datoveho souboru
% 
% Autori: Mgr. & Mgr. Petr Dluhos, RNDr. Eva Janousova
% Masarykova univerzita, Brno
% 05/09/2014
%
% ----------------------------------------------------------------------- %

%nacteme datovy soubor (dokument Microsoft excel)
[num, txt, raw] = xlsread('Data_neuro.xlsx',1);

%num - matice s numerickymi hodnotami
%txt - pole s textovymi hodnotami
%raw - pole se vsemi daty ve forme textu

%vybereme data popisujici spravne zarezeni subjektu do trid
tridy = num(:, 3);

%a 2 priznaky, podle kterych budeme klasifikovat - objem hipokampu a
%amygdaly
data = num(:, 23:24);


% % nebo alternativne 6 priznaku: pohlavi (0/1), MMSE, objem hipokampu a
% %amygdaly, thalamu a pallida
% data = num(:, [9 15 23 24 25 26]); - nevykresli se obrazky, protoze to
% jde jen pro 2D dataset

%vymazeme tridu 2 a ponechame jen skupiny 1 a 3 (kontorly a AD)
data(tridy == 2, :) = [];
tridy(tridy == 2) = [];

%zjistime celkovy pocet subjektu
pocetSubjektu = length(tridy)

%rozdelime mnozinu vsech subjektu na trenovaci a testovaci v pomeru zhruba
%80:20 a ulozime si je do zvlastnich promennych, aby v budoucnu nedoslo k
%zamene.
holdout = cvpartition(pocetSubjektu,'Holdout',0.2)
dataTrenovaciHO = data(holdout.training, :);
dataValidacniHO = data(holdout.test, :);

tridyTrenovaciHO = tridy(holdout.training);
tridyValidacniHO = tridy(holdout.test);

%10-fold krizova validace - rozdelime mnozinu na 10 skupin, pak na 9 ucime
%a na 1 testujeme
cv = cvpartition(pocetSubjektu,'kfold',10)

for i = 1:10
    dataTrenovaci{i} = data(cv.training(i), :);
    dataValidacni{i} = data(cv.test(i), :);

    tridyTrenovaci{i} = tridy(cv.training(i));
    tridyValidacni{i} = tridy(cv.test(i)); 
end    


%% Metoda podpurnych vektoru (SVM - support vector machine)

%nezavisla mnozina

%natrenujeme klasifikator a vizualizujeme si situaci pomoci obrazku.
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true')


%alternativne lze menit hodnotu parametru C (boxconstraint), vychozi je C=1:
C = 0.1;
figure
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true', 'boxconstraint', C)

C = 1;
figure
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true', 'boxconstraint', C)

C = 10;
figure
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true', 'boxconstraint', C)



%pokud metoda nekonverguje (casto pro vysoke hodnoty C), lze zvetsit pocet iteraci SMO:
C = 1000;
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true', 'boxconstraint', C)


options = statset('Display', 'iter', 'MaxIter', 200000)
%a pote pouzit strukturu options na svtupu svmtrain
SVMStruct = svmtrain(dataTrenovaciHO,tridyTrenovaciHO,'showplot','true', 'boxconstraint', C, 'options', options)

klasifikaceTrenovaci = svmclassify(SVMStruct, dataTrenovaciHO);
spatneTrenovaci = sum(klasifikaceTrenovaci ~= tridyTrenovaciHO)
chybaTrenovaci = spatneTrenovaci / length(tridyTrenovaciHO)

klasifikaceValidacni = svmclassify(SVMStruct,dataValidacniHO,'showplot','true');
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO)
chybaValidacni = spatneValidacni / length(tridyValidacniHO)


%krizova validace
%postupne ucime na 9 a testujeme na 1
for i = 1:10
        
    figure
    SVMStruct = svmtrain(dataTrenovaci{i},tridyTrenovaci{i},'showplot','true');
    
    klasifikaceTrenovaci = svmclassify(SVMStruct, dataTrenovaci{i});
    spatneTrenovaci = sum(klasifikaceTrenovaci ~= tridyTrenovaci{i});
    chybyTrenovaci(i) = spatneTrenovaci / length(tridyTrenovaci{i});

    
    klasifikaceValidacni = svmclassify(SVMStruct,dataValidacni{i},'showplot','true');
    spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacni{i});
    chybyValidacni(i) = spatneValidacni / length(tridyValidacni{i});

end 

chybyTrenovaci
chybyValidacni

%prumerne chyby v 10 klasifikacich na trenovaci a validacni mnozine
mean(chybyTrenovaci)
mean(chybyValidacni)



%% Bayesuv klasifikator

[klasifikaceValidacni, err, POSTERIOR, logp, coeff] = classify(dataValidacniHO,dataTrenovaciHO,tridyTrenovaciHO,'quadratic');

%vizualizace
figure
hold on
gscatter(dataTrenovaciHO(:, 1),dataTrenovaciHO(:, 2),tridyTrenovaciHO,'rb','o',5,'off')
K = coeff(1,2).const;
L = coeff(1,2).linear;
Q = coeff(1,2).quadratic;
% Function to compute K + L*v + v'*Q*v for multiple vectors
% v=[x;y]. Accepts x and y as scalars or column vectors.
f = @(x,y) K + [x y]*L + sum(([x y]*Q) .* [x y], 2);
h2 = ezplot(f,[5000 8000 2000 3500]);
set(h2,'Color','m','LineWidth',2)
axis([5000 8000 2000 3500])
xlabel('Objem hipokampu')
ylabel('Objem amygdaly')
title('{\bf Klasifikace pacientu s AD a zdravych kontrol}')
legend('kontroly', 'pacienti')

gscatter(dataValidacniHO(:, 1),dataValidacniHO(:, 2),tridyValidacniHO,'yg','.',15,'off')

% chyba klasifikace
err

%posteriorni pravdepodobnosti zarazeni do jednotlivych trid (kazdy subjekt
%je jeden sloupec
POSTERIOR'


%krizova validace
%postupne ucime na 9 a testujeme na 1
for i = 1:10
        
    [klasifikaceValidacni, err, POSTERIOR, logp, coeff] = classify(dataValidacni{i},dataTrenovaci{i},tridyTrenovaci{i},'quadratic');

    
    chybyValidacni(i) = err;

end 


chybyValidacni

%prumerna chyba v 10 klasifikacich
mean(chybyValidacni)


%% Klasifikace podle minimalni vzdalenosti

%Euklidovska
for i = 1:15
  model = fitcknn(dataTrenovaciHO, tridyTrenovaciHO, 'NumNeighbors', i, 'NSMethod', 'exhaustive', 'Distance', 'euclidean');

  klasifikaceValidacni = model.predict(dataValidacniHO);
  spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
  chybaValidacni = spatneValidacni / length(tridyValidacniHO);
  
  fprintf('KNN, Euklidova metrika, n=%d - chyba na validacni mnozine: %3.4f\n', i, chybaValidacni);
end

fprintf('\n');


%Manhattanska
for i = 1:15
  model = fitcknn(dataTrenovaciHO, tridyTrenovaciHO, 'NumNeighbors', i, 'NSMethod', 'exhaustive', 'Distance', 'cityblock');

  klasifikaceValidacni = model.predict(dataValidacniHO);
  spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
  chybaValidacni = spatneValidacni / length(tridyValidacniHO);
  
  fprintf('KNN, manhattanska metrika, n=%d - chyba na validacni mnozine: %3.4f\n', i, chybaValidacni);
end

fprintf('\n');

%Cebysevova
for i = 1:15
  model = fitcknn(dataTrenovaciHO, tridyTrenovaciHO, 'NumNeighbors', i, 'NSMethod', 'exhaustive', 'Distance', 'chebychev');

  klasifikaceValidacni = model.predict(dataValidacniHO);
  spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
  chybaValidacni = spatneValidacni / length(tridyValidacniHO);
  
  fprintf('KNN, Cebysevova metrika, n=%d - chyba na validacni mnozine: %3.4f\n', i, chybaValidacni);
end

fprintf('\n');

%Canberrska
%cm = @(XI, XJ) (sqrt(sum(bsxfun(@minus, XJ, XI).^2, 2))); %euklidovska
cm = @(XI, XJ) (sum( abs(bsxfun(@minus, XJ, XI)) ./ bsxfun(@plus, abs(XJ), abs(XI)), 2)); %suma (|xi - yi| / (|xi| + |xj|))
%nebude vsak fungovat zcela spravne, pokud citatel i jmenovatel bude 0  


for i = 1:15
  model = fitcknn(dataTrenovaciHO, tridyTrenovaciHO, 'NumNeighbors', i, 'NSMethod', 'exhaustive', 'Distance', cm);

  klasifikaceValidacni = model.predict(dataValidacniHO);
  spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
  chybaValidacni = spatneValidacni / length(tridyValidacniHO);
  
  fprintf('KNN, Canberrska metrika, n=%d - chyba na validacni mnozine: %3.4f\n', i, chybaValidacni);
end


%% average linkage = metoda prumerne vazby

%Euklidovska metrika

%vypocitame vzdalenosti mezi kazdym bodem z validacni a kazdym bodem z trenovaci mnoziny
D = pdist2(dataValidacniHO, dataTrenovaciHO, 'euclidean');

%zobrazime matici
figure
imagesc(D)

%vypocitame prumernou vzdalenost k jednotlivym mnozinam
meanDist1 = mean(D(:, tridyTrenovaciHO == 1), 2);
meanDist3 = mean(D(:, tridyTrenovaciHO == 3), 2);

%vizualizace prumernych vzdalenosti vsech bodu z validacni mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(meanDist1 > meanDist3) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);

fprintf('\n metoda prumerne vazby, Euklidovska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);




%Manhattanska metrika

%vypocitame vzdalenosti mezi kazdym bodem z validacni a kazdym bodem z trenovaci mnoziny
D = pdist2(dataValidacniHO, dataTrenovaciHO, 'cityblock');

%vypocitame prumernou vzdalenost k jednotlivym mnozinam
meanDist1 = mean(D(:, tridyTrenovaciHO == 1), 2);
meanDist3 = mean(D(:, tridyTrenovaciHO == 3), 2);

%vizualizace prumernych vzdalenosti vsech bodu z validacni mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(meanDist1 > meanDist3) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);

fprintf('\n metoda prumerne vazby, manhattanska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);



%Cebysevova metrika

%vypocitame vzdalenosti mezi kazdym bodem z validacni a kazdym bodem z trenovaci mnoziny
D = pdist2(dataValidacniHO, dataTrenovaciHO, 'chebychev');

%vypocitame prumernou vzdalenost k jednotlivym mnozinam
meanDist1 = mean(D(:, tridyTrenovaciHO == 1), 2);
meanDist3 = mean(D(:, tridyTrenovaciHO == 3), 2);

%vizualizace prumernych vzdalenosti vsech bodu z validacni mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(meanDist1 > meanDist3) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n metoda prumerne vazby, Cebysevova metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);



%Canberrska metrika

%vypocitame vzdalenosti mezi kazdym bodem z validacni a kazdym bodem z trenovaci mnoziny
D = pdist2(dataValidacniHO, dataTrenovaciHO, cm);

%vypocitame prumernou vzdalenost k jednotlivym mnozinam
meanDist1 = mean(D(:, tridyTrenovaciHO == 1), 2);
meanDist3 = mean(D(:, tridyTrenovaciHO == 3), 2);

%vizualizace prumernych vzdalenosti vsech bodu z validacni mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(meanDist1 > meanDist3) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n metoda prumerne vazby, Canberrska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);



%% centroidova metoda

%vypocteme centroidy
C1 = mean(dataTrenovaciHO(tridyTrenovaciHO == 1, :))
C3 = mean(dataTrenovaciHO(tridyTrenovaciHO == 3, :))

figure
hold on
gscatter(dataTrenovaciHO(:, 1),dataTrenovaciHO(:, 2),tridyTrenovaciHO,'rb','o',5,'off')
xlabel('Objem hipokampu')
ylabel('Objem amygdaly')
title('{\bf Centoidu jednotlivych skupin}')
legend('kontroly', 'pacienti')
plot(C1(1), C1(2), 'r.', 'MarkerSize',30);
plot(C3(1), C3(2), 'b.', 'MarkerSize',30);


%Euklidovska metrika

%vypocteme vzdalenosti mezi centroidy a vsemi body z validacni mnoziny
D = pdist2([C1; C3], dataValidacniHO,'euclidean');

figure
hold on
plot(D(1, :), 'r'); %vzdalenost bodu k centroidu kontrol
plot(D(2, :), 'b'); %vzdalenost bodu k centroidu AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(D(1, :) > D(2, :)) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n centroidova metoda, Euklidovska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);

%Manhattanska metrika metrika

%vypocteme vzdalenosti mezi centroidy a vsemi body z validacni mnoziny
D = pdist2([C1; C3], dataValidacniHO,'cityblock');

figure
hold on
plot(D(1, :), 'r'); %vzdalenost bodu k centroidu kontrol
plot(D(2, :), 'b'); %vzdalenost bodu k centroidu AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(D(1, :) > D(2, :)) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n centroidova metoda, Manhattanska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);



%Cebysevova metrika

%vypocteme vzdalenosti mezi centroidy a vsemi body z validacni mnoziny
D = pdist2([C1; C3], dataValidacniHO,'chebychev');

figure
hold on
plot(D(1, :), 'r'); %vzdalenost bodu k centroidu kontrol
plot(D(2, :), 'b'); %vzdalenost bodu k centroidu AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(D(1, :) > D(2, :)) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n centroidova metoda, Cebysevova metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);




%Canberrska metrika

%vypocteme vzdalenosti mezi centroidy a vsemi body z validacni mnoziny
D = pdist2([C1; C3], dataValidacniHO, cm);

figure
hold on
plot(D(1, :), 'r'); %vzdalenost bodu k centroidu kontrol
plot(D(2, :), 'b'); %vzdalenost bodu k centroidu AD

klasifikaceValidacni = ones(length(tridyValidacniHO), 1);
klasifikaceValidacni(D(1, :) > D(2, :)) = 3;
spatneValidacni = sum(klasifikaceValidacni ~= tridyValidacniHO);
chybaValidacni = spatneValidacni / length(tridyValidacniHO);
fprintf('\n centroidova metoda, Canberrska metrika\n - chyba na validacni mnozine: %3.4f\n', chybaValidacni);


