Test script updated

test/eeg_import_test.py

@@ -0,0 +1,163 @@
+# Get Paths
+from glob import glob
+
+# EEG package
+from mne import pick_types
+from mne.io import read_raw_edf
+
+import os
+import numpy as np
+
+#%%
+"""
+# Get file paths
+PATH = '/Users/jimmy/data/PhysioNet' #'/rigel/pimri/users/xh2170/data2/data' #PATH = './data/'
+SUBS = glob(os.path.join(PATH, 'S[0-9]*'))
+FNAMES = sorted([x[-4:] for x in SUBS])
+"""
+
+PATH = '/Users/jimmy/data/PhysioNet' #'/rigel/pimri/users/xh2170/data2/data' #PATH = './data/' '..\data\BCI2000'
+SUBS = glob(os.path.join(PATH, 'S[0-9]*'))
+FNAMES = sorted([x[-4:] for x in SUBS])
+#filePath = "..\data\BCI2000\S001" #E:\hxf\Faculty\ColumbiaUniversity\research\myLab\RAs\EEG\Euiyoung(Jimmy)Chung\project\EEG\Demo\data\BCI2000\S001
+#dataPath = filePath + "\S001R04.edf"
+
+try:
+    FNAMES.remove('S089')
+except:
+    pass
+
+#input: a list of subject name, window time in seconds
+#output: shape(X)=[sample#,channel#, windowSize], sample# depends on the window size and step size (defined inside the funciton below), e.g., we have 100 time points, windowSize=10, stepSize=5, we will have 19 samples
+def get_data_forTesting(subj_num, count, interval, epoch_sec=0.0625):
+    """ Import from edf files data and targets in the shape of 3D tensor
+    
+        Output shape: (Trial*Channel*TimeFrames)
+        
+        Some edf+ files recorded at low sampling rate, 128Hz, are excluded. 
+        Majority was sampled at 160Hz.
+        
+        epoch_sec: time interval for one segment of mashes
+        """
+
+    # Event codes mean different actions for two groups of runs
+    run_type_0 = '02'.split(',')
+    run_type_1 = '04,08,12'.split(',')
+    run_type_2 = '06,10,14'.split(',')
+
+    # Initiate X, y
+    X = []
+    y = []
+
+
+    # fixed numbers
+    nChan = 64 
+    sfreq = 160
+    sliding = epoch_sec/2 #step size
+    timeFromCue = 0.5 #exclude 0.5 seconds data right after cue, subjective value here
+
+    # Sub-function to assign X and X, y
+    def append_X(n_segments, old_x):
+        '''This function generate a tensor for X and append it to the existing X'''
+        def window(n):
+            windowStart = int(timeFromCue*sfreq) + int(sfreq*sliding*n)
+            windowEnd = int(timeFromCue*sfreq) + int(sfreq*sliding*(n+2))
+            return [windowStart, windowEnd]
+        new_x = old_x + [data[:, window(n)[0]: window(n)[1]] for n in range(n_segments)\
+                 if data[:, window(n)[0]:window(n)[1]].shape==(nChan, int(sfreq*epoch_sec))]
+        return new_x
+
+    def append_X_Y(run_type, event, old_x, old_y):
+        '''This function seperate the type of events 
+        (refer to the data descriptitons for the list of the types)
+        Then assign X and Y according to the event types'''
+        # Number of sliding windows
+        n_segments = int(event[1]/epoch_sec) #event[1] is the total runnting time for a single run : 4.1s in most cases
+
+        # Instantiate new_x, new_y
+        new_y = old_y
+        new_x = old_x
+
+        # y assignment
+        if run_type == 1:
+            if event[2] == 'T1':
+                new_y = old_y + [1]*n_segments
+                new_x = append_X(n_segments, old_x)
+
+            elif event[2] == 'T2':
+                new_y = old_y + [2]*n_segments
+                new_x = append_X(n_segments, old_x)
+
+        if run_type == 2:
+            if event[2] == 'T1':
+                new_y = old_y + [3]*n_segments
+                new_x = append_X(n_segments, old_x)
+
+            elif event[2] == 'T2':
+                new_y = old_y + [4]*n_segments
+                new_x = append_X(n_segments, old_x)
+
+        return new_x, new_y
+
+    # Iterate over subj_num: S001, S002, S003...
+    for i, subj in enumerate(subj_num):
+        # Return completion rate
+        # Get file names
+        fnames = glob(os.path.join(PATH, subj, subj+'R*.edf'))
+        fnames = [name for name in fnames if name[-6:-4] in run_type_0+run_type_1+run_type_2]
+
+    for i, fname in enumerate(fnames):
+
+            # Import data into MNE raw object
+            raw = read_raw_edf(fname, preload=True, verbose=False)
+
+            picks = pick_types(raw.info, eeg=True)
+
+            if raw.info['sfreq'] != 160:
+                print('{} is sampled at 128Hz so will be excluded.'.format(subj))
+                break
+
+            # High-pass filtering
+            raw.filter(l_freq=1, h_freq=None, picks=picks)
+
+            # Get annotation
+            try:
+                events = raw.find_edf_events()
+            except:
+                continue
+            # Get data
+            data = raw.get_data(picks=picks)
+
+            # Number of this run
+            which_run = fname[-6:-4]
+
+            """ Assignment Starts """ 
+            # run 1 - baseline (eye closed)
+            if which_run in run_type_0:
+
+                # Number of sliding windows
+                n_segments = int((raw.n_times/(epoch_sec*sfreq)))
+
+                # Append 0`s based on number of windows
+                y.extend([0]*n_segments)
+                X = append_X(n_segments, X)
+
+            # run 4,8,12 - imagine opening and closing left or right fist    
+            elif which_run in run_type_1:
+
+                for i, event in enumerate(events):
+                    X, y = append_X_Y(run_type=1, event=event, old_x=X, old_y=y)
+
+            # run 6,10,14 - imagine opening and closing both fists or both feet
+            elif which_run in run_type_2:
+
+                for i, event in enumerate(events):         
+                    X, y = append_X_Y(run_type=2, event=event, old_x=X, old_y=y)
+
+    X = np.stack(X)
+    y = np.array(y)
+
+    how_many_slice = int(interval//epoch_sec)*2 - 1
+    start_slice = int((2*interval//epoch_sec-1)*count)
+
+    return X[start_slice:start_slice+how_many_slice,:,:], y[start_slice:start_slice+how_many_slice]

test/eeg_preprocessing_test.py

@@ -0,0 +1,70 @@
+import numpy as np
+from sklearn.preprocessing import scale
+
+#%%
+def convert_mesh(X):
+
+    mesh = np.zeros((X.shape[0], X.shape[2], 10, 11, 1))
+    X = np.swapaxes(X, 1, 2)
+
+    # 1st line
+    mesh[:, :, 0, 4:7, 0] = X[:,:,21:24]
+
+    # 2nd line
+    mesh[:, :, 1, 3:8, 0] = X[:,:,24:29]
+
+    # 3rd line
+    mesh[:, :, 2, 1:10, 0] = X[:,:,29:38]
+
+    # 4th line
+    mesh[:, :, 3, 1:10, 0] = np.concatenate((X[:,:,38].reshape(-1, X.shape[1], 1),\
+                                          X[:,:,0:7], X[:,:,39].reshape(-1, X.shape[1], 1)), axis=2)
+
+    # 5th line
+    mesh[:, :, 4, 0:11, 0] = np.concatenate((X[:,:,(42, 40)],\
+                                        X[:,:,7:14], X[:,:,(41, 43)]), axis=2)
+
+    # 6th line
+    mesh[:, :, 5, 1:10, 0] = np.concatenate((X[:,:,44].reshape(-1, X.shape[1], 1),\
+                                        X[:,:,14:21], X[:,:,45].reshape(-1, X.shape[1], 1)), axis=2)
+
+
+    # 7th line
+    mesh[:, :, 6, 1:10, 0] = X[:,:,46:55] 
+
+    # 8th line
+    mesh[:, :, 7, 3:8, 0] = X[:,:,55:60] 
+
+    # 9th line
+    mesh[:, :, 8, 4:7, 0] = X[:,:,60:63] 
+
+    # 10th line
+    mesh[:, :, 9, 5, 0] = X[:,:,63] 
+
+    return mesh
+
+#input:...JImported data by eeg_import.py - dim (Trial*Channel*windowSize)
+#output:..JSamples as 2D meshes - dim (Trial*windowSize*height*width*1) -> the last 1 is purely to apply CNN
+
+#%%
+#input X, shape(X)=[sample#,channel#, windowSize], sample# depends on the window size and step size, e.g., we have 100 time points, windowSize=10, stepSize=5, we will have 19 samples
+#step1: one hot encoding for lables
+#step2: shuffle data X
+#step3: standarize X, Z score
+#step4: convert X to 2D matrix, i.e., 10x11 size matrix
+#return: shape(X)= [sample,windowSize,matrixHeight,matrixWidth], label y, shape(y)=[sample#,classes#]
+def prepare_data_forTesting(X, return_mesh = True):
+
+    # Z-score Normalization
+    def scale_data(X):
+        shape = X.shape
+        for i in range(shape[0]):
+            X[i,:, :] = scale(X[i,:, :])
+        return X
+
+    if return_mesh:
+        X_preprocessed  = convert_mesh(scale_data(X))
+
+    return X_preprocessed
+
+

test/eeg_testLive16.py

@@ -0,0 +1,156 @@
+'''
+Author: Euiyoung Chung
+Data: Oct. 10, 2018
+Descriptions:
+
+This script is writted to test the model trained by 'eeg_main.py' for new real time eeg data.
+Device: Emotiv Epoc +
+Number of channel: 14
+Sampling rate: 128 or 256
+Frequency range: 0.16 - 43Hz
+
+'''
+
+# Load in libraries
+from mne.io import read_raw_edf
+from mne import pick_types
+
+from keras.models import load_model 
+from sklearn.preprocessing import scale
+
+import numpy as np
+import time
+
+import signal
+import sys
+
+from collections import defaultdict
+
+# Path setting
+filePath = "./data"
+dataPath = filePath + "/test.edf"
+modelPath = "./model"
+
+# EEG info
+sfreq = 128 #or 256
+
+# (1)Load in data
+def loadData(interval, dataPath = dataPath):
+    raw = read_raw_edf(dataPath, preload=True, verbose=False)
+    eegPicks = pick_types(raw.info, eeg=True)
+    raw.filter(l_freq=1, h_freq=None, picks=eegPicks)
+    return raw.get_data(picks=eegPicks)[-sfreq*interval:,:]
+
+# (2)Preprocess data - 2D meshes + normalization + sliding windows 
+# Assume the input array shape [None, 14]
+# Reshape the input into [None, 10, 6, 6,1] - [batch, window, width, height, 1]
+
+def preprocess(data):
+
+    def scaler(data):
+        ''' The data should be shaped as
+        [None, 14] 
+        '''
+        return scale(data, axis=1)
+
+    def slidingWindow(data):
+
+        n_window = data.shape[0]//10 - 1
+        n_frame = 10
+        n_channel = 14
+
+        newData = np.zeros((n_window, n_frame, n_channel))
+
+        for i in range(n_window):
+            newData[i, :, :] = data[5*i:5*i+10, :]
+
+        return newData
+
+    def mesh(data):
+        n_window = data.shape[0]
+        n_frame = 10
+        height = width = 6
+
+        newData = np.zeros((n_window, n_frame, height, width))
+
+        newData[:, :, 0, np.array([1, 4])] =
+        newData[:, :, 1, np.array([0,2,3,5])] =
+        newData[:, :, 2, np.array([1, 4])] =
+        newData[:, :, 4, np.array([0, 5])] =
+        newData[:, :, 5, np.array([1, 4])] = 
+
+        return newData
+
+    return mesh(slidingWindow(scaler(data)))
+
+# (3)Load in model
+def loadModel(path = modelPath + "/model.h5"):
+    model = load_model(path)
+    return model
+
+model = loadModel()
+
+# (4)Predict - assuming we know the true values
+'''
+if __name__ == "__main__":
+	y_pred = []	
+	for i in testData.shape[0]:
+	 	y_pred.append(model.predict(test).argmax(axis = 1))
+
+	score = classification_report(y_test, y_pred, digits=4, output_dict=True)
+	print(score)
+'''
+
+# (5) Predict - realtime
+total_running = int(sys.argv[2]) #seconds
+interval = int(sys.argv[1]) #seconds
+
+class MyTimer(object):
+    """
+    https://pythonadventures.wordpress.com/2012/12/08/terminate-a-script-after-x-seconds/
+    Similar to Qt's QTimer. Call a function in a specified time.
+ 
+    Time is given in sec. Usage:
+    mt = MyTimer()
+    mt.singleShot(<sec>, <function_name>)
+ 
+    After setting it, you can still disable it:
+    mt.disable()
+ 
+    If you call it several times, any previously scheduled alarm
+    will be canceled (only one alarm can be scheduled at any time).
+    """
+    def singleShot(self, sec, func):
+        self.f = func
+        signal.signal(signal.SIGALRM, self.handler)
+        signal.alarm(sec)
+
+    def handler(self, *args):
+        self.f()
+
+    def disable(self):
+        signal.alarm(0)
+
+def printPred(model, interval):
+
+    while True:
+        data = loadData(interval)
+        data = preprocess(data)
+
+        print(model.predict(data).argmax(axis=1))
+
+        time.sleep(interval)
+
+
+if __name__ == "__main__":
+    mt = MyTimer().singleShot(total_running, sys.exit)        
+    printPred(model, interval)
+
+
+
+
+
+
+
+
+