% Skript s ukazkami klasifikace obrazovych dat mozku; validace
% provedena pomoci LOOCV (leave-one-out cross-validation)
%
% by Eva Koritakova, 01/12/2016, modifikovano 21/11/2017 a 30/10/2018.


%% 0. NACTENI A PREDZPRACOVANI DAT.

% Nastaveni cesty k souborum:
cd 'c:\Users\koritakova\Documents\Vyuka\Vicerozmerky a AKD\Analyza a klasifikace dat\2018-podzim\AKD_predn-07_cviceni\Data_klasifikace\'

% Nacteni souboru se jmeny obrazovych souboru subjektu a zarazenim do skupin:
[num,text]=xlsread('gm_identifikatory.xls','List1'); % num...ciselne hodnoty, text...textove hodnoty
names = text(2:end,1); % jmena souboru s obrazy jednotlivych subjektu
group = num(:,1); % zarazeni subjektu do skupin
inxP = find(group==1); % identifikatory pacientu
inxC = find(group==0); % identifikatory kontrol
clear num text % odstraneni nepotrebnych promennych
group(group==0)=2; % prekodovani vektoru identifikatoru skupin, protoze logisticka regrese vyzaduje nenulove hodnoty

% Nacteni masky sede hmoty mozku (GM) - pro odstraneni mimomozkovych a dalsich sktruktur:
maska=load('gm_mask_2D.mat'); % nacteni masky
maska=maska.gm_mask_2D; % vyber matice s maskou ze strukturovaneho pole
zobraz(maska) % zobrazeni masky GM

% Nacteni 2-D obrazu prvniho pacienta na ukazku:
i=1;
img=load([cell2mat(names(inxP(i),1)),'_rez50.mat']); % nacteni obrazu
zobraz(img.rez) % zobrazeni 2-D GM obrazu prvniho pacienta

% Nacteni vsech obrazovych dat (2-D rezy obrazu mozku):
X=zeros(length(inxP)+length(inxC),sum(sum(maska==1))); % alokace datove matice, kde v radcich budou odmaskovane obrazy jednotlivych subjektu
for i=1:length(inxP) % pro pacienty:
    img=load([cell2mat(names(inxP(i),1)),'_rez50.mat']); % nacteni obrazu pacienta
    X(i,:)=img.rez(maska==1); % odmaskovani obrazu a ulozeni do radku matice X
end
for i=1:length(inxC) % pro kontroly:
    img=load([cell2mat(names(inxC(i),1)),'_rez50.mat']); % nacteni obrazu kontrolniho subjektu
    X(i+length(inxP),:)=img.rez(maska==1); % odmaskovani obrazu a ulozeni do radku matice X
end



%% REDUKCE DAT POMOCI MANNOVA-WHITNEYOVA TESTU + RUZNE KLASIFIKATORY

group_klasif1=zeros(length(group),5); % alokace vektoru s indexy skupin urcenymi klasifikatorem
sel=zeros(length(group),size(X,2)); % alokace vektoru s vybranymi pixely v kazde iteraci
for i=1:length(group)
    i
    x_test=X(i,:); % testovaci subjekt
    X_train=X(~ismember(1:length(group),i),:);
    group_train=group(~ismember(1:length(group),i));
    pvals=zeros(1,size(X_train,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
    for j=1:size(X_train,2) % pro vsechny pixely:
        pvals(j)=ranksum(X_train(group_train==1,j),X_train(group_train==2,j)); % vypocet p-hodnoty pomoci Mannova-Whitneyova testu - srovnani hodnot pacientu a kontrol v danem pixelu
    end
    [pID,pN] = FDR(pvals,0.05); % False Discovery Rate (FDR) korekce
    [sum(pvals<0.05),sum(pvals<pID)] % vypis statisticky vyznamnych pixelu pred a po FDR korekci
    sel(i,:)=double(pvals<pID);
    X_train_red=X_train(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci
    x_test_red=x_test(pvals<pID); % redukce dat testovaciho subjektu
    
    % KLASIFIKATORY:
    SVMStruct = fitcsvm(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif1(i,1)=predict(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru   
    %SVMStruct = svmtrain(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu) - ve starsich verzich Matlabu
    %group_klasif1(i,1)=svmclassify(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru - ve starsich verzich Matlabu
    
    [group_klasif1(i,2),~,~,~,~]=classify(x_test_red,X_train_red,group_train,'diaglinear'); % naivni Bayesuv klasifikator
    
    B = TreeBagger(50,X_train_red,group_train); % nahodne lesy (random forests) - natrenovani
    group_klasif1(i,3) = str2num(cell2mat(predict(B,x_test_red))); % nahodne lesy (random forests) - klasifikace testovaciho subjektu
    
    [B] = mnrfit(X_train_red,group_train); % logistická regrese - natrenovani
    pihat = mnrval(B,x_test_red); % logistická regrese - klasifikace testovaciho subjektu
    if pihat(2)>=0.5
        group_klasif1(i,4)=2;
    else
        group_klasif1(i,4)=1;
    end
    
    model = fitcknn(X_train_red,group_train,'NumNeighbors',5,'NSMethod','exhaustive','Distance','euclidean'); % k-NN s 5 sousedy a Euklidovskou vzdalenosti - natrenovani
    group_klasif1(i,5) = model.predict(x_test_red); % k-NN s 5 sousedy a Euklidovskou vzdalenosti - klasifikace testovaciho subjektu
    
end

% Vyhodnoceni uspesnosti klasifikace:
effic=zeros(size(group_klasif1,2),7); % alokace matice s mirami uspesnosti
for i=1:size(group_klasif1,2)
    eq =(group_klasif1(:,i)==group); % a logicky vector: 0..klasifikace je nespravna, 1..klasifikace je spravna

    TP=sum(eq(group==1));    % skutecne pozitivni vysledky (true positive)
    FN=sum(group==1)-TP;     % falesne negativni vysledky (false negative)
    TN=sum(eq(~(group==1))); % skutecne negativni vysledky (true negative)
    FP=sum(~(group==1))-TN;  % falesne pozitivni vysledky (false positive)

    accur = (TP+TN) / (TP+TN+FP+FN); % spravnost (accuracy)
    sens = TP / (TP+FN);             % senzitivita (sensitivity)
    spec = TN / (TN+FP);             % specificita (specificity)

    effic(i,:)=[TP FN TN FP accur sens spec]; % poskladani do jednoho vektoru
end
% effic =
%    19.0000    5.0000   17.0000    7.0000    0.7500    0.7917    0.7083 - SVM
%    18.0000    6.0000   19.0000    5.0000    0.7708    0.7500    0.7917 - naivni Bayesuv klasifikator
%    18.0000    6.0000   18.0000    6.0000    0.7500    0.7500    0.7500 - klasifikacni stromy
%    15.0000    9.0000   16.0000    8.0000    0.6458    0.6250    0.6667 - logisticka regrese
%    18.0000    6.0000   17.0000    7.0000    0.7292    0.7500    0.7083 - k-NN s 5 sousedy a Euklidovskou vzdalenosti

% Srovnani uspesnosti klasifikace k-NN s nahodnou klasifikaci (pomoci jednovyberoveho binomickeho testu):
if sum(group==1)>sum(group==2)
    PA0=sum(group==1)/length(group); % uspesnost nahodne klasifikace, pokud je pocet subjektu v 1. skupine vetsi nez ve 2. skupine
else
    PA0=sum(group==2)/length(group); % uspesnost nahodne klasifikace, pokud je pocet subjektu ve 2. skupine vetsi nebo roven 1. skupine
end
p_val=myBinomTest(effic(5,1)+effic(5,3),length(group),PA0,'two'); % prvni prvek je TN+TP pro k-NN klasifikator
p_val % 0.0021 - tzn. uspesnost klasifikace statisticky vyznamne lepsi nez nahodna klasifikace

% Srovnani uspesnoti klasifikace mezi k-NN a naivnim Bayesovym klasifikatorem (pomoci McNemarova testu):
eq2=(group_klasif1(:,[2,5])==(group*ones(1,2))); % a logic matrix: 0..classification is incorrect, 1..classification is correct
N01=sum((eq2(:,1)==0).*(eq2(:,2)==1)); % calculating N01
N10=sum((eq2(:,1)==1).*(eq2(:,2)==0)); % calculating N10
if N01+N10~=0
    chi2=((abs(N01-N10)-1)^2)/(N01+N10); % testova statistika
    p=1-chi2cdf(chi2,1); % p-hodnota
else
    p=1; % p-hodnota
end
p % 0.4795 - tzn. neni statisticky vyznamny rozdil v uspesnosti mezi k-NN a naivnim Bayesovym klasifikatorem

% Vizualizace vyznamnych pixelu - pixely, ktere byly statisticky vyznamne
% nejmene v polovine iteraci:
vyzn=maska;
vyzn(maska==1)=1+((sum(sel)/size(sel,1))>=0.5);
zobraz(vyzn)
