% Skript s ukazkami redukce a klasifikace obrazovych dat mozku se zamerenim
% na ruzne metody validace.
%
% by Eva Koritakova (Janousova), 11/09/2014, modifikovano 29/11/2016.


%% 0. 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. REKDUKCE A KLASIFIKACE DAT - VALIDACE POMOCI RESUBSTITUCE.
% (tzn. trenovaci a testovaci soubor je stejny)

% redukce dat vyberem statisticky vyznamnych pixelu pomoci Mannova-Whitneyova testu:
pvals=zeros(1,size(X,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
for i=1:size(X,2) % pro vsechny pixely:
    pvals(i)=ranksum(X(group==1,i),X(group==0,i)); % vypocet p-hodnoty pomoci Mannova-Whitneyova testu - srovnani hodnot pacientu a kontrol v danem pixelu
end
sum(pvals<0.05) % pocet statisticky vyznammych pixelu
[pID,pN] = FDR(pvals,0.05); % False Discovery Rate (FDR) korekce
sum(pvals<pID) % pocet statisticky vyznammych pixelu po FDR korekci
X_red=X(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci

% zobrazeni statisticky vyznamych pixelu:
vyzn=maska;
vyzn(maska==1)=1+(pvals<pID);
zobraz(vyzn)

% klasifikace redukovanych obrazovych dat:
SVMStruct = svmtrain(X_red,group); % natrenovani SVM klasifikatoru
%SVMStruct = svmtrain(X_red,group,'boxconstraint', C); %alternativne lze menit hodnotu parametru C (boxconstraint), vychozi je C=1:
group_klasif = svmclassify(SVMStruct, X_red); % klasifikace dat pomoci SVM klasifikatoru
group_klasif=group_klasif(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 100%; 100%; 100%



%% 2. REDUKCE A KLASIFIKACE DAT - VALIDACE POMOCI METODY LEAVE-ONE-OUT (NESPRAVNA VARIANTA)
% (nespravna varianta proto, ze krosvalidace je delana jen pro klasifikaci
% a ne i pro redukci dat)

% redukce dat vyberem statisticky vyznamnych pixelu pomoci Mannova-Whitneyova testu:
pvals=zeros(1,size(X,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
for i=1:size(X,2) % pro vsechny pixely:
    pvals(i)=ranksum(X(group==1,i),X(group==0,i)); % vypocet p-hodnoty pomoci Mannova-Whitneyova testu - srovnani hodnot pacientu a kontrol v danem pixelu
end
sum(pvals<0.05) % pocet statisticky vyznammych pixelu
[pID,pN] = FDR(pvals,0.05); % False Discovery Rate (FDR) korekce
sum(pvals<pID) % pocet statisticky vyznammych pixelu po FDR korekci
X_red=X(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci

% klasifikace redukovanych obrazovych dat:
group_klasif2=zeros(1,length(group)); % alokace vektoru s vysledky klasifikace (tzn. indexy skupin urcenymi klasifikatorem)
for i=1:length(group)
    i
    x_test_red=X_red(i,:); % vyber testovaciho subjektu ze zredukovanych dat
    X_train_red=X_red(~ismember(1:length(group),i),:); % vyber trenovacich subjektu ze zredukovanych dat
    group_train=group(~ismember(1:length(group),i)); % vyber identifikatoru trenovacich 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
% 81.3%; 79.2%; 83.3%



%% 3. REDUKCE A KLASIFIKACE DAT + VALIDACE POMOCI METODY LEAVE-ONE-OUT (SPRAVNA VARIANTA)

group_klasif3=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_klasif3(i)=svmclassify(SVMStruct, x_test_red); % 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
% 72.3%; 70.8%; 75.0%

% POZNAMKA: je patrne z vypisu statisticky vyznamnych pixelu, ze je pocet
% statisticky vyznamnych pixelu v kazde iteraci hodne variabilni, a tudiz
% se musi delat vyber pixelu jen pomoci trenovacich dat, jinak dostaneme
% zkreslene vysledky!

% 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:
vyzn3=maska;
vyzn3(maska==1)=1+((sum(sel)/size(sel,1))>=0.5);
zobraz(vyzn3)



%% 4. REDUKCE A KLASIFIKACE DAT + VALIDACE POMOCI 10-NASOBNE KRIZOVE VALIDACE

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
    SVMStruct = svmtrain(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif4(cv.test(i))=svmclassify(SVMStruct, X_test_red); % klasifikace dat pomoci SVM klasifikatoru
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%; 70.8%; 79.2%
% 66.7%; 66.7%; 66.7%
% 75.0%; 66.7%; 83.3%

% 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:
vyzn4=maska;
vyzn4(maska==1)=1+((sum(sel)/10)>=0.5); % ve jmenovateli je 10, protoze bylo 10 iteraci
zobraz(vyzn4)



%% 5. REDUKCE A KLASIFIKACE DAT + HOLD-OUT

holdout = cvpartition(length(group),'Holdout',0.25); % rozdeleni subjektu na 75% trenovacich a 25% testovacich

X_test=X(holdout.test,:); % testovaci subjekty
X_train=X(holdout.training,:); % trenovaci subjekty
group_train=group(holdout.training);
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
if isempty(pID)
    display('Zadne pixely nejsou vyznamne po FDR korekci, klasifikaci tedy nelze provest')
else
    [sum(pvals<0.05),sum(pvals<pID)] % vypis statisticky vyznamnych pixelu pred a po FDR korekci
    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_klasif5=svmclassify(SVMStruct, X_test_red); % klasifikace dat pomoci SVM klasifikatoru
    group_test=group(holdout.test);
    spravnost = sum(group_klasif5 == group_test) / length(group_test) % vypocet spravnosti klasifikace
    senzitivita = sum(group_klasif5(group_test==1) == group_test(group_test==1)) / sum(group_test==1) % vypocet senzitivity klasifikace
    specificita = sum(group_klasif5(group_test==0) == group_test(group_test==0)) / sum(group_test==0) % vypocet specificity klasifikace
end
% pri kazdem spusteni dostavame trochu jine vysledky podle toho, jak se
% data nahodne rozdeli na trenovaci a testovaci; vysledky napr.:
% 83.3%; 75.0%; 100.0%
% 50.0%; 50.0%; 50.0%
% 75.0%; 80.0%; 71.4%
% 66.7%; 62.5%; 75.0%
% nebo dokonce v nekterych pripadech ani klasifikace neprobehne, protoze po
% FDR korekci neni vyznamny ani jeden pixel

% zobrazeni statisticky vyznamych pixelu:
if isempty(pID)
    display('Zadne pixely nejsou vyznamne po FDR korekci, klasifikaci tedy nelze provest')
else
    vyzn5=maska;
    vyzn5(maska==1)=1+(pvals<pID);
    zobraz(vyzn5)
end
