main_classify.py

# -*- coding: utf-8 -*-
"""
Evaluating five different adaptive decomposition methods for EEG signal seizure detection and classification
Vinícius R. Carvalho,  Márcio F.D. Moraes, Antônio P. Braga, Eduardo M.A.M. Mendes
Programa de Pós-Graduação em Engenharia Elétrica – Universidade Federal de Minas Gerais – Av. Antônio Carlos 6627, 31270-901, Belo Horizonte, MG, Brasil.

This script first reads the files generated by main_feats.py, with features extracted from each decomposed method
by EMD, EEMD, CEEMDAN, EWT or VMD. Then, it splits training/test samples according to k-folds method, trains and evaluates
each dataset with different classifiers. This is done 10x, resulting in performance tables with mean and std values.


@author: Vinicius Carvalho
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time
import seaborn as sns
from sklearn.preprocessing import StandardScaler, MultiLabelBinarizer, LabelBinarizer
from sklearn.decomposition import PCA
from sklearn.neural_network import MLPClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import roc_auc_score,accuracy_score, confusion_matrix
from sklearn.feature_selection import SelectKBest,chi2, RFE
from skfeature.function.similarity_based import fisher_score 


def fun_classify(inputFile, groupsSel, FeatSelect, Nfeats,scaleFeats = 1):
    """
    AllStatsMean, AllStatsSTD = fun_classify(inputFile, groupsSel, FeatSelect, Nfeats)
    inputFile: the .csv file containt feature tables
    groups: The selected groups to classify. Full set is ["S","F","Z","N","O"],
    but ["S","F","Z"] are of most interest for the article (ictal, inter-ictal and normal EEG)
    FeatSelect: feature selection method: PCA, RFE, fisher or none
    Nfeats: number of selected features
    Returns:
    AllStatsMean: mean performance values
    AllStatsSTD: standard deviation of performance values  
    """
    #reads input features
    dfFeats = pd.read_csv(inputFile, sep=',',header=0)

    #only selected groups
    dfFeats = dfFeats[dfFeats["Group"].isin(groupsSel)]
    if "decTaime" in dfFeats:
        x = dfFeats.iloc[:, 2:]#ignores decomposition method execution time
    else:
        x = dfFeats.iloc[:, 1:]
    y = dfFeats.iloc[:, 0].values
    if scaleFeats:#scale feats?
        x = StandardScaler().fit_transform(x)
    #Feature selection
    if x.shape[1] > Nfeats:
        #RFE
        if FeatSelect == "RFE":
            rfeModel = SVC(kernel="linear", C=0.025,probability = True,gamma = 'scale')
            rfeSelect = RFE(rfeModel,n_features_to_select = Nfeats)
            rfe_fit = rfeSelect.fit(x, y)
            x = x[:,rfe_fit.support_]
            
        if FeatSelect == "PCA":
            pca = PCA(n_components=Nfeats)
            x = pca.fit_transform(x)
        
        if FeatSelect == "fisher":
            fisherScore = fisher_score.fisher_score(x, y)
            idx = fisher_score.feature_ranking(fisherScore)
            x = x[:,idx[:Nfeats]]

    names = ["KNN", "Linear SVM", "RBF SVM", "GPC", "MLP"] 
    
    classifiers = [
        KNeighborsClassifier(3),
        SVC(kernel="linear", C=0.025,probability = True,gamma = 'scale'),
        SVC(probability = True,gamma = 'scale'),
        GaussianProcessClassifier(1.0 * RBF(1.0)),
        MLPClassifier(alpha=1,max_iter = 200)]

    #initialize performance variable
    AllStats = {}
    AllStatsMean = {}   
    AllStatsSTD = {}   
    
    for name in names:
        AllStats[name] = {"Accuracy":np.zeros([realizations,K_folds]),
            "SensitivityMean":np.zeros([realizations,K_folds]),
            "SpecificityMean":np.zeros([realizations,K_folds]),
            "AUC_Mean":np.zeros([realizations,K_folds]),
            "SensitivityIctal":np.zeros([realizations,K_folds]),
            "SpecificityIctal":np.zeros([realizations,K_folds]),
            "AUC_Ictal":np.zeros([realizations,K_folds]),
            "TTtimes":np.zeros([realizations,K_folds])}  
        AllStatsMean[name] = {"Accuracy":0.,"SensitivityMean":0.,
                    "SpecificityMean":0,"AUC_Mean":0.,"SensitivityIctal":0.,
                    "SpecificityIctal":0.,"AUC_Ictal":0.,"TTtimes":0.}
        AllStatsSTD[name] = {"Accuracy":0.,"SensitivityMean":0.,
                    "SpecificityMean":0,"AUC_Mean":0.,"SensitivityIctal":0.,
                    "SpecificityIctal":0.,"AUC_Ictal":0., "TTtimes":0.}    
        #for each realization
    for i in range(realizations):
        skf = StratifiedKFold(n_splits=K_folds,shuffle = True) #5-fold validation
        
        for tupTemp,ki in zip(skf.split(x, y),range(K_folds)):
            train_idx, test_idx = tupTemp[0],tupTemp[1]
            X_train, X_test = x[train_idx], x[test_idx]
            y_train, y_test = y[train_idx], y[test_idx]   
            for name, clf in zip(names, classifiers):   #for each classifier
                tic = time.time()#check training/testing time of each classifier
                #Fit model and predict
                modelFit = clf.fit(X_train, y_train)
                yPredicted = modelFit.predict(X_test)
                probsTest = modelFit.predict_proba(X_test)
                toc = time.time()
                # AUC -  #ictal class as positive 
                if len(np.unique(y)) > 2:
                    AUCs = roc_auc_score(LabelBinarizer().fit_transform(y_test), probsTest, average = None)
                else:  
                    AUCs = roc_auc_score(y_test, probsTest[:,1], average = None)
                #Sensitivity and Specificity
                cMatrix = confusion_matrix(y_test, yPredicted)        
                FP = cMatrix.sum(axis=0) - np.diag(cMatrix)  
                FN = cMatrix.sum(axis=1) - np.diag(cMatrix)
                TP = np.diag(cMatrix)
                TN = cMatrix.sum() - (FP + FN + TP)
                # Sensitivity
                TPR = TP/(TP+FN)
                # Specificity or true negative rate
                TNR = TN/(TN+FP) 
                #fill performance variable
                AllStats[name]["Accuracy"][i,ki] = accuracy_score(y_test, yPredicted)
                AllStats[name]["SensitivityMean"][i,ki] = np.mean(TPR)
                AllStats[name]["SpecificityMean"][i,ki] = np.mean(TNR)
                AllStats[name]["SensitivityIctal"][i,ki] = TPR[0]
                AllStats[name]["SpecificityIctal"][i,ki] = TNR[0]                
                AllStats[name]["AUC_Mean"][i,ki] = np.mean(AUCs)
                AllStats[name]["TTtimes"][i,ki] = toc-tic
                if len(np.unique(y)) > 2:
                    AllStats[name]["AUC_Ictal"][i,ki] = AUCs[0]
    AllStatsDF = [0]*len(names)
    for idx, name in enumerate(names):    
        for istat in AllStats[name].keys():
            AllStats[name][istat] = np.mean(AllStats[name][istat],axis = 1)
            AllStatsMean[name][istat] = np.mean(AllStats[name][istat])
            AllStatsSTD[name][istat] = np.std(AllStats[name][istat])  
        AllStatsDF[idx] = pd.DataFrame.from_dict(AllStats[name])
        AllStatsDF[idx]["Nmodes"] = Nmodes
        AllStatsDF[idx]["Classifier"] = name
        
    return pd.DataFrame.from_dict(AllStatsMean),pd.DataFrame.from_dict(AllStatsSTD),  pd.concat(AllStatsDF)

#%% main script
    
dBase = "NSC_ND" #which dataset (NSC-ND or BonnDataset)

groups = {"BonnDataset":["S","F","Z"],
          "NSC_ND":["ictal","interictal","preictal"]} #groups to include
    
#Define parameters
Nmodes_list = [2,3,4,5,6,7,8] #number of modes for decomposition [EMD, EWT, VMD] - check for input files
numFeatures = 20 #if FeatSelect = "PCA" or "RFE", chooses the number of used features

realizations = 10#number of realizations
K_folds = 10 #number of k-folds
scaleFeats = 1 #if 1, z-scores all features
FeatSelect = "RFE" #Feature selection method: PCA, kbest (TODO), RFE, fisher. if 0, uses all features

RESULTS_Export = [0]*len(Nmodes_list)
RESULTS_DF = [0]*len(Nmodes_list)
Nmodes = 1
RESULTS_Orig_mean,RESULTS_Orig_std, RESULTS_Orig_DF = fun_classify('%s/ORIGFeatsWelch.csv'%dBase, groups[dBase], FeatSelect, 11)
RESULTS_Orig_DF["Method"] = "Orig"

for idxN, Nmodes in enumerate(Nmodes_list):
    RESULTS = {"EMD":0,"EEMD":0,"CEEMDAN":0,"EWT":0,"VMD":0}
    RESULTS_std = {"EMD":0,"EEMD":0,"CEEMDAN":0,"EWT":0,"VMD":0}
    RESULTS["EMD"],RESULTS_std["EMD"],DF1 = fun_classify('%s/EMDFeatsWelch_%dModes.csv'%(dBase,Nmodes), groups[dBase], FeatSelect, numFeatures,scaleFeats)
    print("EMD ok")
    RESULTS["EEMD"],RESULTS_std["EEMD"],DF2 = fun_classify('%s/EEMDFeatsWelch_%dModes.csv'%(dBase,Nmodes), groups[dBase], FeatSelect, numFeatures,scaleFeats)
    print("EEMD ok")
    RESULTS["CEEMDAN"],RESULTS_std["CEEMDAN"],DF3 = fun_classify('%s/CEEMDANFeatsWelch_%dModes.csv'%(dBase,Nmodes), groups[dBase], FeatSelect, numFeatures,scaleFeats)
    print("CEEMDAN ok")
    RESULTS["EWT"],RESULTS_std["EWT"],DF4 = fun_classify('%s/EWTFeatsWelch_%dModes.csv'%(dBase,Nmodes), groups[dBase], FeatSelect, numFeatures,scaleFeats)
    print("EWT ok")
    RESULTS["VMD"],RESULTS_std["VMD"],DF5 = fun_classify('%s/VMDFeatsWelch_%dModes.csv'%(dBase,Nmodes), groups[dBase], FeatSelect, numFeatures,scaleFeats)
    print("VMD ok")
    DF1["Method"] = "EMD"
    DF2["Method"] = "EEMD"
    DF3["Method"] = "CEEMDAN"
    DF4["Method"] = "EWT"
    DF5["Method"] = "VMD"
    RESULTS_DF[idxN] = pd.concat([DF1,DF2,DF3,DF4,DF5], ignore_index=True)
    RESULTS_Export[idxN] = {"EMD":0,"EEMD":0,"CEEMDAN":0,"EWT":0,"VMD":0,"Orig":0}
    for ii in list(RESULTS.keys()):
        RESULTS[ii].loc["AUC_Mean"] /= 100
        AA = RESULTS[ii].applymap(lambda x:'{:s}'.format("%.2f"%(round(100*x,2))))
        BB = RESULTS_std[ii].applymap(lambda x:'{:s}'.format("%.2f"%(round(100*x,2))))
        RESULTS_Export[idxN][ii] = AA + '+-' + BB
        
    
    RESULTS_Export[idxN]["ALL"] = pd.DataFrame.append(RESULTS_Export[idxN]["EMD"],RESULTS_Export[idxN]["EEMD"])
    RESULTS_Export[idxN]["ALL"] = pd.DataFrame.append(RESULTS_Export[idxN]["ALL"],RESULTS_Export[idxN]["CEEMDAN"])
    RESULTS_Export[idxN]["ALL"] = pd.DataFrame.append(RESULTS_Export[idxN]["ALL"],RESULTS_Export[idxN]["EWT"])
    RESULTS_Export[idxN]["ALL"] = pd.DataFrame.append(RESULTS_Export[idxN]["ALL"],RESULTS_Export[idxN]["VMD"])
    #RESULTS_Export[idxN]["ALL"] = pd.DataFrame.append(RESULTS_Export[idxN]["ALL"],RESULTS_Export[idxN]["Orig"])

RESULTS_DF = pd.concat(RESULTS_DF,ignore_index = True)
RESULTS_DF = pd.concat([RESULTS_Orig_DF, RESULTS_DF],ignore_index = True)   

RESULTS_DF.to_pickle("{}/{}{}_ClassResults{}.pkl".format(dBase,FeatSelect,numFeatures,Nmodes_list))
RESULTS_DF["Accuracy"] *= 100

AA = RESULTS_Orig_mean.applymap(lambda x:'{:s}'.format("%.2f"%(round(100*x,2))))
BB = RESULTS_Orig_std.applymap(lambda x:'{:s}'.format("%.2f"%(round(100*x,2))))
RESULTS_Export[-1]["Orig"] = AA + '+-' + BB

#plot accuracy as a function of the number of modes for a single classifier (RBF/Linear SVM, KNN, MLP, GPC)
figClassifier = "RBF SVM"
fig, ax1 = plt.subplots(figsize=(3.543, 2.2))
sns.barplot(ax = ax1, x = "Nmodes",y = "Accuracy",
            data = RESULTS_DF[RESULTS_DF["Classifier"] == figClassifier], 
            hue = "Method",
            errwidth = 1)
plt.title("%s - %s"%(figClassifier,dBase))
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)  
plt.xlabel("Number of modes")
if dBase == "NSC_ND":
    plt.ylim([55,100])#adjust y axis
else:
    plt.ylim([85,100])
plt.ylabel("ACC (%)")
ax1.locator_params(axis='y', nbins=5)
plt.tight_layout()
sns.plotting_context("paper")
ax1.legend().set_visible(False)
plt.savefig('%s%d_%s_ACCxNmodes_%s.pdf'%(FeatSelect,numFeatures,figClassifier,dBase), dpi = 300)
plt.savefig('%s%d_%s_ACCxNmodes_%s.png'%(FeatSelect,numFeatures,figClassifier,dBase), dpi = 300)