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


%% NACTENI A PREDZPRACOVANI DAT.

% Nastaveni cesty k souborum:
cd 'c:\Users\koritakova\Documents\Vyuka\DSAN02\2018 jaro\DSAN02_predn08_klasifikace2\prakticka_cast_1\Data_a_pomocne_soubory\'

% 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



%% 1. REDUKCE DAT POMOCI MANNOVA-WHITNEYOVA TESTU

group_klasif1=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:length(group)
    i
    x_test=X(i,:); % vyber testovaciho subjektu
    X_train=X(~ismember(1:length(group),i),:); % vyber trenovacich subjektu
    group_train=group(~ismember(1:length(group),i)); % vyber identifikatoru trenovacich subjektu
    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); % vytvoreni matice s vybranymi vyznamnymi pixely (slouzi pro vykresleni vyzn. pixelu)
    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
    SVMStruct = svmtrain(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif1(i)=svmclassify(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru
end
group_klasif1=group_klasif1(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif1 == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif1(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif1(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 72.3%; 70.8%; 75.0%

% vizualizace vyznamnych pixelu - v kazde iteraci dostaneme jiny pocet
% vyznamnych pixelu; pri vizualizaci muzeme napr. vykreslit ty pixely, 
% ktere byly statisticky vyznamne nejmene v polovine iteraci:
vyzn1=maska;
vyzn1(maska==1)=1+((sum(sel)/size(sel,1))>=0.5);
zobraz(vyzn1)



%% 2. REDUKCE DAT POMOCI PCA

group_klasif2=zeros(1,length(group)); % alokace vektoru s indexy skupin urcenymi klasifikatorem
for i=1:length(group)
    i
    X_train=X(~ismember(1:length(group),i),:); % vyber trenovacich subjektu
    group_train=group(~ismember(1:length(group),i)); % vyber identifikatoru trenovacich subjektu
    
    [coeff,~,~] = pca(X_train); % coeff - matice vlastnich vektoru (vl. vektory jsou ve sloupcich)
    X_train_red=X_train*coeff; % vynasobeni dat matici vlastnich vektoru
    X_train_mean=mean(X_train_red); % vypocet prumeru z kazde hlavni komponenty
    X_train_red=X_train_red-ones(size(X_train_red,1),1)*X_train_mean; % centrovani (tzn. odecteni prumeru z kazde hlavni komponenty)
    
    x_test=X(i,:); % vyber testovaciho subjektu
    x_test_red=x_test*coeff; % redukce dat testovaciho subjektu
    x_test_red=x_test_red-ones(size(x_test_red,1),1)*X_train_mean; % centrovani dat testovaciho subjektu
    
    SVMStruct = svmtrain(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif2(i)=svmclassify(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru
end
group_klasif2=group_klasif2(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif2 == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif2(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif2(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 66.7%; 66.7%; 66.7%


%% 3. REDUKCE DAT POMOCI PCA + MANNOVA-WHITNEYOVA TESTU

group_klasif3=zeros(1,length(group)); % alokace vektoru s indexy skupin urcenymi klasifikatorem
for i=1:length(group)
    i
    X_train=X(~ismember(1:length(group),i),:); % vyber trenovacich subjektu
    group_train=group(~ismember(1:length(group),i)); % vyber identifikatoru trenovacich subjektu
    
    % PCA:
    [coeff,~,~] = pca(X_train); % coeff - matice vlastnich vektoru (vl. vektory jsou ve sloupcich)
    X_train_red=X_train*coeff; % vynasobeni dat matici vlastnich vektoru
    X_train_mean=mean(X_train_red); % vypocet prumeru z kazde hlavni komponenty
    X_train_red=X_train_red-ones(size(X_train_red,1),1)*X_train_mean; % centrovani (tzn. odecteni prumeru z kazde hlavni komponenty)
    
    % M-W test:
    pvals=zeros(1,size(X_train_red,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
    for j=1:size(X_train_red,2) % pro vsechny pixely:
        pvals(j)=ranksum(X_train_red(group_train==1,j),X_train_red(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 - prilis prisna, neni tam nic vyznamneho, proto pID nastavime na 0.05
    pID=0.05;
    X_train_red2=X_train_red(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci
    
    % redukce testovacich dat:
    x_test=X(i,:); % vyber testovaciho subjektu
    x_test_red=x_test*coeff; % redukce dat testovaciho subjektu
    x_test_red=x_test_red-ones(size(x_test_red,1),1)*X_train_mean; % centrovani dat testovaciho subjektu
    x_test_red2=x_test_red(pvals<pID); % redukce dat testovaciho subjektu
    
    % klasifikace:
    SVMStruct = svmtrain(X_train_red2,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif3(i)=svmclassify(SVMStruct, x_test_red2); % klasifikace dat pomoci SVM klasifikatoru
end
group_klasif3=group_klasif3(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif3 == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif3(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif3(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 64.6%; 75.0%; 54.2% - u techto dat vyber vyzn. komponent pomohl jen ke
% zvyseni senzitivity (nikoliv ke zvyseni celkove spravnosti ani
% specificity)

