% UKOL - vypocitejte spravnost, senzitivitu a specificitu klasifikace 
%        pomoci k-NN klasifikatoru s k=9 s Hammingovou (manhattanskou)  
%        metrikou vzdalenosti a 10-fold krosvalidaci
%      - redukci provedte selekci voxelu statisticky vyznamnych dle
%        Mannova-Whitneyova testu po FDR korekci
%      - napoveda: vyuzijte skripty 02_Prehled_metod_validace.m a 
%        03_Prehled_metod_klasifikace.m


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

% 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



%% 2. REDUKCE DAT POMOCI MANNOVA-WHITNEYOVA TESTU + k-NN KLASIFIKATOR + 10-NASOBNA KRIZOVA VALIDACE

%rng(10) % zafixovani nahodneho generatoru, abychom pri kazdem spusteni dostali stejny vysledek
cv = cvpartition(length(group),'kfold',10); % rozdeleni subjektu do 10 skupin
group_klasif4=zeros(1,length(group)); % 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:cv.NumTestSets
    i
    X_test=X(cv.test(i),:); % testovaci subjekty
    X_train=X(cv.training(i),:); % trenovaci subjekty
    group_train=group(cv.training(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==0,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
    model = fitcknn(X_train_red,group_train,'NumNeighbors',9,'NSMethod','exhaustive','Distance','cityblock'); % k-NN s 9 sousedy a Hammingovou (manhattanskou) vzdalenosti - natrenovani
    group_klasif4(cv.test(i)) = model.predict(X_test_red); % k-NN s 9 sousedy a Hammingovou (manhattanskou) vzdalenosti - klasifikace testovaciho subjektu
end
group_klasif4=group_klasif4(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif4 == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif4(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif4(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% pri kazdem spusteni dostavame trochu jine vysledky podle toho, jak se
% data nahodne rozdeli do 10 skupin; vysledky napr.:
% 75.0%; 75.0%; 75.0%
% 77.1%; 75.0%; 79.2%
% 75.0%; 83.3%; 66.7%
% 75.0%; 79.2%; 70.8%
% 68.8%; 75.0%; 62.5%
