%Klasifikace pomocí dat z velkého datového souboru
%by Petr Dluhos 2014

%nacteme datovy soubor (je potrebne mit nastavenou cestu do adresare)
[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 zarazeni subjektu do trid
tridy = num(:, 3);

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


% % 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 (kontroly 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 (hold-out metoda) a ulozime si je do zvlastnich promennych, aby v 
%budoucnu nedoslo k zamene
holdout = cvpartition(pocetSubjektu,'Holdout',0.2)
dataTrenovaciHO = data(holdout.training, :);
dataTestovaciHO = data(holdout.test, :);
tridyTrenovaciHO = tridy(holdout.training);
tridyTestovaciHO = tridy(holdout.test);

%10-fold krizova validace - rozdelime mnozinu na 10 skupin, pak vzdy na 9 z 
%nich ucime a na 1 testujeme
cv = cvpartition(pocetSubjektu,'kfold',10)
for i = 1:10
    dataTrenovaci{i} = data(cv.training(i), :);
    dataTestovaci{i} = data(cv.test(i), :);

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


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

% 1. hold-out metoda:

%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)


%klasifikace trenovacich dat (dostaneme nadhodnocene vysledky)
klasifikaceTrenovaci = svmclassify(SVMStruct, dataTrenovaciHO);
spatneTrenovaci = sum(klasifikaceTrenovaci ~= tridyTrenovaciHO)
chybaTrenovaci = spatneTrenovaci / length(tridyTrenovaciHO)

%klasifikace testovacich dat
klasifikaceTestovaci = svmclassify(SVMStruct,dataTestovaciHO,'showplot','true');
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO)
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO)
%spravnost muzeme vypocitat jako 1-chybaTestovaci

% 2. krizova validace:
%postupne ucime na 9 a testujeme na 1 podsouboru
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});

    klasifikaceTestovaci = svmclassify(SVMStruct,dataTestovaci{i},'showplot','true');
    spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaci{i});
    chybyTestovaci(i) = spatneTestovaci / length(tridyTestovaci{i});

end 

chybyTrenovaci
chybyTestovaci

%prumerne chyby v 10 klasifikacich na trenovaci a Testovaci mnozine
mean(chybyTrenovaci)
mean(chybyTestovaci)



%% Bayesuv klasifikator

% 1. hold-out metoda:

[klasifikaceTestovaci, err, POSTERIOR, logp, coeff] = classify(dataTestovaciHO,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(dataTestovaciHO(:, 1),dataTestovaciHO(:, 2),tridyTestovaciHO,'yg','.',15,'off')

% chyba klasifikace
err

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


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

end 

chybyTestovaci

%prumerna chyba v 10 klasifikacich
mean(chybyTestovaci)



%% Klasifikace podle minimalni vzdalenosti

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

  klasifikaceTestovaci = model.predict(dataTestovaciHO);
  spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
  chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
  
  fprintf('KNN, Euklidova metrika, n=%d - chyba na Testovaci mnozine: %3.4f\n', i, chybaTestovaci);
end
fprintf('\n');

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

  klasifikaceTestovaci = model.predict(dataTestovaciHO);
  spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
  chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
  
  fprintf('KNN, manhattanska metrika, n=%d - chyba na Testovaci mnozine: %3.4f\n', i, chybaTestovaci);
end
fprintf('\n');

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

  klasifikaceTestovaci = model.predict(dataTestovaciHO);
  spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
  chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
  
  fprintf('KNN, Cebysevova metrika, n=%d - chyba na Testovaci mnozine: %3.4f\n', i, chybaTestovaci);
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|)) ; neni ale osetren pripad, kdy citatel i jmenovatel jsou 0  
for i = 1:15
  model = fitcknn(dataTrenovaciHO, tridyTrenovaciHO, 'NumNeighbors', i, 'NSMethod', 'exhaustive', 'Distance', cm);

  klasifikaceTestovaci = model.predict(dataTestovaciHO);
  spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
  chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
  
  fprintf('KNN, Canberrska metrika, n=%d - chyba na Testovaci mnozine: %3.4f\n', i, chybaTestovaci);
end


%% average linkage = metoda prumerne vazby

%Euklidovska metrika

%vypocitame vzdalenosti mezi kazdym bodem z Testovaci a kazdym bodem z trenovaci mnoziny
D = pdist2(dataTestovaciHO, 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 Testovaci mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(meanDist1 > meanDist3) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n metoda prumerne vazby, Euklidovska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Manhattanska metrika

%vypocitame vzdalenosti mezi kazdym bodem z Testovaci a kazdym bodem z trenovaci mnoziny
D = pdist2(dataTestovaciHO, 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 Testovaci mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(meanDist1 > meanDist3) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n metoda prumerne vazby, manhattanska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Cebysevova metrika

%vypocitame vzdalenosti mezi kazdym bodem z Testovaci a kazdym bodem z trenovaci mnoziny
D = pdist2(dataTestovaciHO, 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 Testovaci mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(meanDist1 > meanDist3) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n metoda prumerne vazby, Cebysevova metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Canberrska metrika

%vypocitame vzdalenosti mezi kazdym bodem z Testovaci a kazdym bodem z trenovaci mnoziny
D = pdist2(dataTestovaciHO, 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 Testovaci mnoziny
figure
hold on
plot(meanDist1, 'r'); %vzdalenost bodu ke kontrolam
plot(meanDist3, 'b'); %vzdalenost bodu k AD

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(meanDist1 > meanDist3) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n metoda prumerne vazby, Canberrska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);



%% 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 Testovaci mnoziny
D = pdist2([C1; C3], dataTestovaciHO,'euclidean');

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

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(D(1, :) > D(2, :)) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n centroidova metoda, Euklidovska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Manhattanska metrika metrika

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

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

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(D(1, :) > D(2, :)) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n centroidova metoda, Manhattanska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Cebysevova metrika

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

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

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(D(1, :) > D(2, :)) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n centroidova metoda, Cebysevova metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);


%Canberrska metrika

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

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

klasifikaceTestovaci = ones(length(tridyTestovaciHO), 1);
klasifikaceTestovaci(D(1, :) > D(2, :)) = 3;
spatneTestovaci = sum(klasifikaceTestovaci ~= tridyTestovaciHO);
chybaTestovaci = spatneTestovaci / length(tridyTestovaciHO);
fprintf('\n centroidova metoda, Canberrska metrika\n - chyba na Testovaci mnozine: %3.4f\n', chybaTestovaci);
