%% Příloha k diplomové práci: ANALÝZA EXPERIMENTÁLNÁCH ZÁZNAMŮ, autor: Maršánová Lucie
%% KLASIFIKACE
%  Skript slouží k realizaci klasifikace EG cyklů do skupin fyziologických cyklů,
%  ischemických cyklů z 5 té minuty ischemie, ischemických cyklů z 10. minuty ischemie 
%  a komorových extrasystol. Při klasifikaci jsou cykly reprezentovány morfologickými
%  parametry (jejich popis viz text DP, jejich výpočet viz další přiložený
%  skript). Cykly je možné reprezentovat několika skupinami parametrů a to
%  s jejich původními hodnotami, nebo normalizovanými. Dále je také možné
%  použít k reprezentaci cyklů parametry vypočtené analýzou hlavních komponent.
%  Klasifikace je realizována 4 klasifikačními metodami a to: metoda
%  nejbližších sousedů, metoda podpůrných vektorů (SVM), diskriminační
%  analýza a Naivní Bayesův klasifikátor. Více viz hlavní text DP.
%%
% 1. část realizuje klasifikaci všemi metodami při všech různých nastaveních
% za použití parametrů získaných analýzou hlavních komponent, 
% 2. část za použití ostatních skupin parametrů - zde je nutné postupně
% jednotlivé skupiny parametrů a konkrétní nastavení klasifikátorů manuálně nastavit
%%
clear all;clc;close all;
load('data.mat') % načtení vstupních dat: param - matice parametrů, paramOut - vektor označující skupiny

% Definice skupin příznaků, které budou použity pro klasifikaci (vybraných
% na základě předchozího statistického testování (Kruskall-Wallis a
% Tukey-Kramer test)
pindAll=[1 2 4 7 15 16 17 21 25 26 27 29 33 39 42 46 47 49 50 51 57 64 65 66 67 69];
pindUnlead=[1 2 4 7];
pindLeadI=[15 16 17 21 25 26 27 29 33 39];
pindLeadII=[42 46 47 49 50 51 57 64 65 66 67 69];

% Definice množin příznků, které použity na testování při 10ti násobné
% křížové validaci - tyto množiny budou tedy vždy pro všechny testované
% metody stejné.
rng(0,'twister'); % nastavení, aby generování začínalo vždy stejně
nfold=10;
cv = cvpartition(paramOut,'kfold',nfold);  
indices = zeros(size(paramOut));
for i=1:nfold
    indices(cv.test(i)) = i;
end

%% 1. Klasifikace při reprezentaci příznaky vytvořenými analýzou hlavních komponent
paramForPCA=param(:,[1:11 13 17:48 50 54:end]); % definice množiny příznaků
accKNNPCA1=[];accKNNPCA5=[];accKNNPCA10=[];accKNNPCA15=[];accKNNPCA20=[];accKNNPCA27=[];
accLDAPCA1=[];accLDAPCA2=[];accLDAPCA3=[];accLDAPCA4=[];accLDAPCA5=[];
accNBPCA1=[]; accNBPCA2=[];
for i=1:nfold % pro první až nfold iteraci křížové validace
    testIdx = (indices == i); trainIdx = ~testIdx; % definice indexů trénovací a testovaí množiny 
    
    ptr=paramForPCA(trainIdx,:); % výběr testovací množiny (dle výše definovanýc indexů)
    mitr=min(ptr); 
    matr=max(ptr); 
    ptrnorm=(ptr-repmat(mitr,size(ptr,1),1))./(repmat(matr,size(ptr,1),1)-repmat(mitr,size(ptr,1),1)); % normalizace treninkové množiny
    
    pte=paramForPCA(testIdx,:); % select testing set
    ptenorm=(pte-repmat(mitr,size(pte,1),1))./(repmat(matr,size(pte,1),1)-repmat(mitr,size(pte,1),1)); % normalizace testovací množiny
    
    [pc,score,latent,tsquare,explained] = princomp(ptrnorm); % výpočet PCA z normalizované treninkové množiny
    
    ptrpca=score(:,1:10); % výběr prvních 10ti hlavních komponent, které budou použity pro další zpracování - klasifikaci
    ptepca=ptenorm*pc(:,1:10); % výpočet prvních 10ti hlavních komponent testovací množiny dat, s využitím vah (load) vypočtených z treninkové mnoziny
    ptepca=ptepca-repmat(mean(ptepca),size(ptepca,1),1); % centrování testovacích PCA
    
    N=sum(testIdx==1); % počet testovacích vektorů
    %% Klasifikace při použití parametrů vypočtených analýou hlavních komponent
    % Metoda nejbližších sousedů, různé nastavení hodnoty k (počet
    % nejbližších sousedů)
    knnClassPCA1 = knnclassify(ptepca, ptrpca, paramOut(trainIdx), 1, 'euclidean');
    knnClassPCA5 = knnclassify(ptepca, ptrpca, paramOut(trainIdx), 5, 'euclidean');
    knnClassPCA10 = knnclassify(ptepca, ptrpca, paramOut(trainIdx),10, 'euclidean');
    knnClassPCA15 = knnclassify(ptepca, ptrpca, paramOut(trainIdx), 15, 'euclidean');
    knnClassPCA20 = knnclassify(ptepca, ptrpca, paramOut(trainIdx), 20, 'euclidean');
    badKNNPCA1 = ~strcmp(knnClassPCA1,paramOut(testIdx)); % výpočet špatné klasifikovaných cyklů (vektorů)
    badKNNPCA5 = ~strcmp(knnClassPCA5,paramOut(testIdx));
    badKNNPCA10 = ~strcmp(knnClassPCA10,paramOut(testIdx));
    badKNNPCA15 = ~strcmp(knnClassPCA15,paramOut(testIdx));
    badKNNPCA20 = ~strcmp(knnClassPCA20,paramOut(testIdx));
    accKNNPCA1(i) = 1-sum(badKNNPCA1)/N; % výpočet úspěšnosti
    accKNNPCA5(i) = 1-sum(badKNNPCA5)/N;
    accKNNPCA10(i) = 1-sum(badKNNPCA10)/N;
    accKNNPCA15(i) = 1-sum(badKNNPCA15)/N;
    accKNNPCA20(i) = 1-sum(badKNNPCA20)/N;
       
    % Diskriminační analýza, různé nastavení diskriminační funkce
    ldaClassPCA1 = classify(ptepca,ptrpca,paramOut(trainIdx),'linear');
    ldaClassPCA2 = classify(ptepca,ptrpca,paramOut(trainIdx),'diaglinear');
    ldaClassPCA3 = classify(ptepca,ptrpca,paramOut(trainIdx),'quadratic');
    ldaClassPCA4 = classify(ptepca,ptrpca,paramOut(trainIdx),'diagquadratic');
    ldaClassPCA5 = classify(ptepca,ptrpca,paramOut(trainIdx),'mahalanobis');
    badLDAPCA1 = ~strcmp(ldaClassPCA1,paramOut(testIdx)); % výpočet špatné klasifikovaných cyklů (vektorů)
    badLDAPCA2 = ~strcmp(ldaClassPCA2,paramOut(testIdx));
    badLDAPCA3 = ~strcmp(ldaClassPCA3,paramOut(testIdx));
    badLDAPCA4 = ~strcmp(ldaClassPCA4,paramOut(testIdx));
    badLDAPCA5 = ~strcmp(ldaClassPCA5,paramOut(testIdx));
    accLDAPCA1(i) = 1-sum(badLDAPCA1)/N; % výpočet úspěšnosti
    accLDAPCA2(i) = 1-sum(badLDAPCA2)/N;
    accLDAPCA3(i) = 1-sum(badLDAPCA3)/N;
    accLDAPCA4(i) = 1-sum(badLDAPCA4)/N;
    accLDAPCA5(i) = 1-sum(badLDAPCA5)/N;
    
    % Naivní Bayesův klasifikátor, různé nastavení typu hustoty
    % pravděpodobnosti
    nbKDPCA1 = NaiveBayes.fit(ptrpca, paramOut(trainIdx)); % trénování klasifikátoru (Gausova hustota pravděpodobnosti)
    nbKDPCA2 = NaiveBayes.fit(ptrpca, paramOut(trainIdx),'dist','kernel'); % trénování klasifikátoru (Kernelova hustota pravděpodobnosti)
    nbClassPCA1 = nbKDPCA1.predict(ptepca); % testování
    nbClassPCA2 = nbKDPCA2.predict(ptepca); % testování      
    badNBPCA1 = ~strcmp(nbClassPCA1,paramOut(testIdx)); % výpočet špatné klasifikovaných cyklů (vektorů)
    badNBPCA2 = ~strcmp(nbClassPCA2,paramOut(testIdx));
    accNBPCA1(i) = 1-sum(badNBPCA1)/N; % výpočet úspěšnosti
    accNBPCA2(i) = 1-sum(badNBPCA2)/N;
end
% SVM klasifikace
% nutne nejprve nainstalovat toolbox viz
% http://www.csie.ntu.edu.tw/~cjlin/libsvm a poté vyndat ze složky
% SVM_pomocne_funkce další pomocné funkce a vložit do stejné, ve které se 
% nachází tento skript, které realizují výpočet za definovaných podmínek 
% křížové validace 
% labels=grp2idx(paramOut);
% opts1 = '-s 0 -t 0 -c 1 -b 1'; % Gausova 
% opts2 = '-s 0 -t 2 -c 1 -b 1'; % jádrová
% accSVMPCA1 = libsvmcrossval_ova_pca(labels, paramForPCA, opts1, nfold);  % funkce, která využívá funkce z toolboxu http://www.csie.ntu.edu.tw/~cjlin/libsvm - je potřeba nainstalovat
% accSVMPCA2 = libsvmcrossval_ova_pca(labels, paramForPCA, opts2, nfold);

% Zobrazení průměrné úspěšnosti (Acc) jednotlivých metod :
accKNNPCA1 = mean(accKNNPCA1)
accKNNPCA5 = mean(accKNNPCA5)
accKNNPCA10 = mean(accKNNPCA10)
accKNNPCA15 = mean(accKNNPCA15)
accKNNPCA20 = mean(accKNNPCA20)
accKNNPCA27 = mean(accKNNPCA27)

accLDAPCA1 = mean(accLDAPCA1)
accLDAPCA2 = mean(accLDAPCA2)
accLDAPCA3 = mean(accLDAPCA3)
accLDAPCA4 = mean(accLDAPCA4)
accLDAPCA5 = mean(accLDAPCA5)

accNBPCA1 = mean(accNBPCA1)
accNBPCA2 = mean(accNBPCA2)

% accSVMPCA1 = mean(accSVMPCA1)
% accSVMPCA2 = mean(accSVMPCA2)

%% 2. Klasifikace cyklů, při reprezentaci 4 různými skupinami příznaků (nutné manuálně nastavit) a s různými nastavením (opět lze manuálně nastavit)

%% metoda nejbližších sousedů 
type='euclidean'; % nastavení typu počítané vzdálenosti
k=1; %1,5,10,15,20 % počet uvažovaných nejbližších sousedů
pind=pindAll;%pindLeadII;%pindLeadI;%pindUnlead; % výběr skupiny příznaků

accInit=zeros(nfold,1); % deficnice proměnné pro ukládání úspěšnosti klasifikace s originálními hodnotami parmametrů
accNorm = zeros(nfold,1); % deficnice proměnné pro ukládání úspěšnosti klasifikace s normalizovanými hodnotami parmametrů

for i=1:nfold % pro každou iteraci
    testIdx = (indices == i); trainIdx = ~testIdx; % select actual testing and training indexes
    
    ptr=param(trainIdx, pind); % select actual training set
    
    mitr=min(ptr); % minimum z treninkoveho setu
    matr=max(ptr); % maximum z treninkoveho setu
    ptrnorm=(ptr-repmat(mitr,size(ptr,1),1))./(repmat(matr,size(ptr,1),1)-repmat(mitr,size(ptr,1),1)); % normalizace treninovací množiny
    
   
    pte=param(testIdx, pind); % výběr aktální testovací množiny
    ptenorm=(pte-repmat(mitr,size(pte,1),1))./(repmat(matr,size(pte,1),1)-repmat(mitr,size(pte,1),1)); % normalizace testovací množiny
       
    
    knnClassInit = knnclassify(pte, ptr, paramOut(trainIdx), k, type); % trénování a klasifikace - originální hodnoty
    knnClassNorm = knnclassify(ptenorm,ptrnorm,paramOut(trainIdx),k,type); % trénování a klasifikace - normalizované hodnoty
        
    badInit = ~strcmp(knnClassInit,paramOut(testIdx)); % výpočet špatné klasifikovaných cyklů (vektorů)
    badNorm = ~strcmp(knnClassNorm,paramOut(testIdx));
    N=sum(testIdx==1);
    accInit(i) = 1-sum(badInit)/N; % výpočet úspěšnosti
    accNorm(i) = 1-sum(badNorm)/N;
        
    [knnResubCMInit,grpOrderInit] = confusionmat(paramOut(testIdx),knnClassInit) % zobrazení matice záměn
    [knnResubCMNorm,grpOrderNorm] = confusionmat(paramOut(testIdx),knnClassNorm)
    end
accKNNInit = mean(accInit) % accuraces
accKNNNorm = mean(accNorm)

%% diskriminační analýza 
type='diagquadratic'; %'diaglinear'; % 'mahalanobis';%'quadratic'; % 'linear';% výběr typu diskriminační funkce
% pind=pindAll;%pindLeadII;%pindLeadI;%pindUnlead; % výběr skupiny příznaků

accInit=zeros(nfold,1);
accNorm = zeros(nfold,1);

for i=1:nfold
    testIdx = (indices == i); trainIdx = ~testIdx;
    
    ptr=param(trainIdx, pind);
    
    mitr=min(ptr); 
    matr=max(ptr);
    meantr=mean(ptr);
    stdtr=std(ptr);
    
    ptrnorm=(ptr-repmat(mitr,size(ptr,1),1))./(repmat(matr,size(ptr,1),1)-repmat(mitr,size(ptr,1),1));
       
    pte=param(testIdx, pind);
    ptenorm=(pte-repmat(mitr,size(pte,1),1))./(repmat(matr,size(pte,1),1)-repmat(mitr,size(pte,1),1));
       
    ldaClassInit = classify(pte,ptr,paramOut(trainIdx),type);
    ldaClassNorm = classify(ptenorm,ptrnorm,paramOut(trainIdx),type);
        
    badInit = ~strcmp(ldaClassInit,paramOut(testIdx));
    badNorm = ~strcmp(ldaClassNorm,paramOut(testIdx));
        
    N=sum(testIdx==1);
    
    accInit(i) = 1-sum(badInit)/N;
    accNorm(i) = 1-sum(badNorm)/N;
        
    [ldaResubCMInit,grpOrderInit] = confusionmat(paramOut(testIdx),ldaClassInit)
    [ldaResubCMNorm,grpOrderNorm] = confusionmat(paramOut(testIdx),ldaClassNorm)
       
end
accLDAInit = mean(accInit)
accLDANorm = mean(accNorm)

%% SVM klasifikace:
% nutne nejprve nainstalovat toolbox viz
% http://www.csie.ntu.edu.tw/~cjlin/libsvm a poté vyndat ze složky
% SVM_pomocne_funkce další pomocné funkce a vložit do stejné, ve které se 
% nachází tento skript, které realizují výpočet za definovaných podmínek 
% křížové validace a normalizace dat

% pind=pindAll;%pindLeadII;%pindLeadI;%pindUnlead; % výběr skupiny příznaků
% labels=grp2idx(paramOut); % popisky
% opts = '-s 0 -t 2 -c 1 -b 1'; % '-s 0 -t 0 -c 1 -b 1'   % výběr jádrové funkce
% 
% data=param(:,pind); % select particular feature set
% acc = libsvmcrossval_ova(labels, data, opts, nfold); 
% fprintf('Acc = %.4f%%\n', 100*mean(acc));
% 
% accn = libsvmcrossval_ova_normal(labels, param, opts, nfold,pind); 
% fprintf('Acc = %.4f%%\n', 100*mean(accn));


%% Naivní Bayesův klasifikátor

% pind=pindAll;%pindLeadII;%pindLeadI;%pindUnlead;% výběr skupiny příznaků
accInit=zeros(nfold,1);
accNorm = zeros(nfold,1);
accStan = zeros(nfold,1);

for i=1:nfold
    testIdx = (indices == i); trainIdx = ~testIdx;
    
    ptr=param(trainIdx, pind);
    
    mitr=min(ptr); 
    matr=max(ptr);
    meantr=mean(ptr);
    stdtr=std(ptr);
    
    ptrnorm=(ptr-repmat(mitr,size(ptr,1),1))./(repmat(matr,size(ptr,1),1)-repmat(mitr,size(ptr,1),1));
    ptrstan=(ptr-repmat(meantr,size(ptr,1),1))./repmat(stdtr,size(ptr,1),1);
    
    pte=param(testIdx, pind);
    ptenorm=(pte-repmat(mitr,size(pte,1),1))./(repmat(matr,size(pte,1),1)-repmat(mitr,size(pte,1),1));
    ptestan=(pte-repmat(meantr,size(pte,1),1))./repmat(stdtr,size(pte,1),1);
    
    % Gausovo rozložení pravděpodobnosti:
%     nbGauInit = NaiveBayes.fit(ptr, paramOut(trainIdx));
%     nbGauNorm = NaiveBayes.fit(ptrnorm, paramOut(trainIdx));
    
%     ldaClassInit = nbGauInit.predict(pte);
%     ldaClassNorm = nbGauNorm.predict(ptenorm);

    
    % Kernelovo rozložení pravděpodobnosti:
    nbKDInit = NaiveBayes.fit(ptr, paramOut(trainIdx),'dist','kernel');
    nbKDNorm = NaiveBayes.fit(ptrnorm, paramOut(trainIdx),'dist','kernel');
    
    
    nbClassInit = nbKDInit.predict(pte);
    nbClassNorm = nbKDNorm.predict(ptenorm);
   
    badInit = ~strcmp(nbClassInit,paramOut(testIdx));
    badNorm = ~strcmp(nbClassNorm,paramOut(testIdx));
  
    N=sum(testIdx==1);
    
    accInit(i) = 1-sum(badInit)/N;
    accNorm(i) = 1-sum(badNorm)/N;
 
    [nbResubCMInit,grpOrderInit] = confusionmat(paramOut(testIdx),nbClassInit);
    [nbResubCMNorm,grpOrderNorm] = confusionmat(paramOut(testIdx),nbClassNorm);
  
end
accNBInit = mean(accInit)
accNBNorm = mean(accNorm)




