% Skript s ukazkami redukce a klasifikace obrazovych dat mozku; validace
% provedena pomoci LOOCV (leave-one-out cross-validation)
%
% by Eva Koritakova, 01/12/2016.


%% NACTENI A PREDZPRACOVANI DAT.

% Nastaveni cesty k souborum:
cd 'c:\Users\koritakova\Documents\Vyuka\Vicerozmerky a AKD\Analyza a klasifikace dat - vyuka\2016-podzim\AKD_predn-10_praktika_Matlab\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

% Nastaveni cesty k funkci recog_effic:
cd 'c:\Users\koritakova\Documents\Vyuka\Vicerozmerky a AKD\Analyza a klasifikace dat - vyuka\2016-podzim\AKD_predn-10_praktika_Matlab\'



%% REDUKCE DAT POMOCI MANNOVA-WHITNEYOVA TESTU

group_klasif1=zeros(length(group),6); % 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 = svmtrain(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif1(i,1)=svmclassify(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru   
    
    [group_klasif1(i,2),~,~,~,~]=classify(x_test_red,X_train_red,group_train,'diaglinear'); % LDA
    
    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 (nelze spocitat kvuli p>n)
    %pihat = mnrval(B,x_test_red); % logistická regrese - klasifikace testovaciho subjektu
    %if pihat(2)>=0.5
    %    group_klasif1(i,4)=2;
    %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
    
    %[group_klasif1(i,6),~,~,~,~] = classify(x_test_red,X_train_red,group_train,'quadratic'); % Bayesuv klasifikator (nelze spocitat kvuli p>n)
end
[effic]=recog_effic(group_klasif1',group,1); % vypocet uspesnosti; posledni prvek udava identifikator pacientu (v techto datech jsou pacienti oznaceni hodnotou 1)
effic

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