{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport scipy.io\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\npreictal_count =0\ninterictal_count=0\nPREICTAL_PATHS = []\nINTERICTAL_PATHS=[]\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        if \"preictal\" in filename and \"Dog\" in filename:\n            preictal_count +=1\n            PREICTAL_PATHS.append(os.path.join(dirname, filename))\n        if \"interictal\" in filename and \"Dog\" in filename:\n            interictal_count +=1\n            INTERICTAL_PATHS.append(os.path.join(dirname, filename))\n            \nprint('  preictal data: ',preictal_count)\nprint('interictal data: ',interictal_count)\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-10T21:29:28.895777Z","iopub.execute_input":"2022-11-10T21:29:28.896218Z","iopub.status.idle":"2022-11-10T21:29:31.593220Z","shell.execute_reply.started":"2022-11-10T21:29:28.896180Z","shell.execute_reply":"2022-11-10T21:29:31.592082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DATA EXPLORATION","metadata":{}},{"cell_type":"markdown","source":"Inspect data stucture","metadata":{}},{"cell_type":"code","source":"mat =scipy.io.loadmat(PREICTAL_PATHS[0])\nprint(mat)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:31.596703Z","iopub.execute_input":"2022-11-10T21:29:31.597584Z","iopub.status.idle":"2022-11-10T21:29:31.833498Z","shell.execute_reply.started":"2022-11-10T21:29:31.597543Z","shell.execute_reply":"2022-11-10T21:29:31.832484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## random plots\nUseful to get an overall idea about the data","metadata":{}},{"cell_type":"code","source":"# Data extraction\nsearch_key = 'segment'\n\n# choose a random ID\ni=np.random.randint(0,len(PREICTAL_PATHS)-1)\npre_mat = scipy.io.loadmat(PREICTAL_PATHS[i])\nsegment = dict(filter(lambda item: search_key in item[0], pre_mat.items()))\n\"\"\"\nIDEA: to extract any segment one should find the element containing the word '_segment_' \ninside the dictionary. To do this you can use a for or adopt a more complex but straightfoward method\n\"\"\"\npreictal_eeg = list(segment.values())[0][0][0][0]\n\nj=np.random.randint(0,len(INTERICTAL_PATHS)-1)\ninter_mat = scipy.io.loadmat(INTERICTAL_PATHS[j])\nsegment = dict(filter(lambda item: search_key in item[0], inter_mat.items()))\ninterictal_eeg = list(segment.values())[0][0][0][0]\n\nprint(f'random preictal file: {PREICTAL_PATHS[i]}')\nprint(f'random interictal file: {INTERICTAL_PATHS[j]}')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:31.834926Z","iopub.execute_input":"2022-11-10T21:29:31.835304Z","iopub.status.idle":"2022-11-10T21:29:32.279040Z","shell.execute_reply.started":"2022-11-10T21:29:31.835262Z","shell.execute_reply":"2022-11-10T21:29:32.278002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Data plotting\nfig0,axs0 = plt.subplots(8,2, figsize=(30,5),sharey=True,sharex=True)\nfig1,axs1 = plt.subplots(8,2, figsize=(30,5),sharey=True,sharex=True)\n\nN_samples = interictal_eeg.shape[1]\nfs = 400\nts = 1/fs\nt0 = np.arange(0,N_samples*ts,ts )\n\nN_samples = preictal_eeg.shape[1]\nt1 = np.arange(0,N_samples*ts,ts )\n\ns=0\nfor r in range(int(interictal_eeg.shape[0]/2)):\n    for c in range(2):\n        axs0[r][c].plot(t0,interictal_eeg[s,:])\n        axs0[r][c].set_title('%ith sensor'%s)\n        axs0[r][c].set_xlim((0,t0[-1]))\n        axs0[r][c].set_xlabel('time (s)')\n        axs0[r][c].set_ylabel('$\\mu Volt$')\n        s+=1\ns=0\n\nfor r in range(int(preictal_eeg.shape[0]/2)):\n    for c in range(2):\n        axs1[r][c].plot(t1,preictal_eeg[s,:])\n        axs1[r][c].set_title('%ith sensor'%s)\n        axs1[r][c].set_xlim((0,t1[-1]))\n        axs1[r][c].set_xlabel('time (s)')\n        axs1[r][c].set_ylabel('$\\mu Volt$')\n        s+=1\n\nfig0.suptitle(f'INTERICTAL EEG ({INTERICTAL_PATHS[j]})')\nfig1.suptitle(f'PREICTAL EEG ({PREICTAL_PATHS[i]})')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:32.280794Z","iopub.execute_input":"2022-11-10T21:29:32.281519Z","iopub.status.idle":"2022-11-10T21:29:37.324274Z","shell.execute_reply.started":"2022-11-10T21:29:32.281479Z","shell.execute_reply":"2022-11-10T21:29:37.323396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del fig0,fig1","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:37.327197Z","iopub.execute_input":"2022-11-10T21:29:37.327889Z","iopub.status.idle":"2022-11-10T21:29:37.334592Z","shell.execute_reply.started":"2022-11-10T21:29:37.327853Z","shell.execute_reply":"2022-11-10T21:29:37.333537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Time frequency analysis","metadata":{}},{"cell_type":"code","source":"!pip install scipy\nfrom scipy import signal","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:37.336277Z","iopub.execute_input":"2022-11-10T21:29:37.336820Z","iopub.status.idle":"2022-11-10T21:29:49.105817Z","shell.execute_reply.started":"2022-11-10T21:29:37.336784Z","shell.execute_reply":"2022-11-10T21:29:49.104627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\nimport sys\nfrom tqdm import tqdm\n\nprint(\"processing preictal data: \")\nZs = 0\nfor i in tqdm(range(len(PREICTAL_PATHS))):\n    pre_mat = scipy.io.loadmat(PREICTAL_PATHS[i])\n    segment = dict(filter(lambda item: search_key in item[0], pre_mat.items()))\n    preictal_eeg = list(segment.values())[0][0][0][0]\n     \n    for channel in range(preictal_eeg.shape[0]): \n        freq,time,Z = scipy.signal.stft(preictal_eeg[channel,:],fs, nperseg=1000)\n        Zs+=np.abs(Z)**2\n    del Z\n    \nZs = Zs/(preictal_eeg.shape[0])\nZs = np.array(Zs)\n    \nplt.pcolormesh(time, freq, Zs, vmin=0)\nplt.ylim((0,30))\nplt.title('preictal spectrum over time ')\nplt.show()\n    \nprint(\"processing interictal data: \")\nfor j in tqdm(range(len(INTERICTAL_PATHS))):\n    j=np.random.randint(0,len(INTERICTAL_PATHS)-1)\n    inter_mat = scipy.io.loadmat(INTERICTAL_PATHS[j])\n    \n    segment = dict(filter(lambda item: search_key in item[0], inter_mat.items()))\n    interictal_eeg = list(segment.values())[0][0][0][0]\n    \n    for channel in range(interictal_eeg.shape[0]): \n        freq,time,Z = scipy.signal.stft(interictal_eeg[channel,:],fs, nperseg=1000)\n        Zs+=np.abs(Z)**2\nZs = Zs/(interictal_eeg.shape[0])\nZs = np.array(Zs)\n\nplt.pcolormesh(time, freq, Zs, vmin=0)\nplt.ylim((0,30))\nplt.title('interictal spectrum over time ')\nplt.show()\n\ndel Z, Zs, freq\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:29:49.111189Z","iopub.execute_input":"2022-11-10T21:29:49.113869Z","iopub.status.idle":"2022-11-10T21:38:36.761377Z","shell.execute_reply.started":"2022-11-10T21:29:49.113826Z","shell.execute_reply":"2022-11-10T21:38:36.760312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems like preictal data have more energy overall however we have less preictal data. This means that higher frequencies are not averaged-off.<br>\nWhat is peculiar are the sudden spikes at low frequences in preictal data. It's still not unique to preictal but it's a sign of what to look for.","metadata":{}},{"cell_type":"code","source":"# if RAM is full it's time to clear some space\nlocal_vars = list(locals().items())\nfor var, obj in local_vars:\n    try:\n        print(var,': ', sys.getsizeof(obj)/1000, 'GB')\n    except:\n        pass","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:38:36.762892Z","iopub.execute_input":"2022-11-10T21:38:36.763237Z","iopub.status.idle":"2022-11-10T21:38:36.771376Z","shell.execute_reply.started":"2022-11-10T21:38:36.763201Z","shell.execute_reply":"2022-11-10T21:38:36.770291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# FRACTAL ANALYSIS\n","metadata":{}},{"cell_type":"markdown","source":"In fractal geometry, the Higuchi dimension (or Higuchi fractal dimension (HFD)) is an approximate value for the box-counting dimension of the graph of a real-valued function or time series.\nFor some types of biological signals, non-linear methods yield qualitatively new information.<br>\n\n## Findings about Fractal Analysis in epilepsy**\n**\"Higuchi Fractal Properties of Onset Epilepsy Electroencephalogram:\"**<br>\nTruong Quang Dang Khoa, Vo Quang Ha, and Vo Van Toi<br>\nhttps://www.hindawi.com/journals/cmmm/2012/461426/ <br>\n\n*[...]the most remarkable aspects of these trends are, during the preictal period, the fractal dimension was relatively high and erratic fluctuates in a small range. However, because Higuchi algorithm is very sensitive to noise [2], especially white noise, the average Fractal dimension of each channel in data is so high and it is so difficult to detect epilepsy. The current difficulty is that we cannot know exactly where the main component from results of ICA is processing. Therefore, we propose using averaging filter for the original data. <br>\nAccording to Figure 6 obtained by the Higuchi algorithm, we can see that fractal coefficient of ECG turned for the worse in the transition from preictal to ictal period. That general pattern was very close to the result of EEG when the fall of fractal figure was marked as the beginning of the seizure. In addition, before the seizure by several minutes, there were two troughs that need to be focus on in anticipation of the seizure.* \n\n\nL. D. lasemidis, J. C. Principe, and J. C. Sackellares, **“Measurement and quantification of Spatio-temporal dynamics of human epileptic seizures,”** in Nonlinear Biomedical Signal Processing: Dynamic Analysis and Modeling, vol. 2 of IEEE Press Series on Biomedical Engineering Metin Akay, Series Editor, 1999.<br>\n\n\n*Perhaps the most exciting discovery to emerge from dynamical analysis of the EEG in \ntemporal lobe epilepsy is that human epileptic seizures are preceded by dynamical \nchanges in the EEG signal. This phenomenon could not be detected by visual inspection \nof the original EEG signal or by other more traditional methods of signal processing. The \nnature of the phenomenon indicates that it may be possible to more accurately localize the \n 26\nepileptogenic focus as well as to predict the onset of a seizure in time to intervene with \nabortive therapy*.\n","metadata":{}},{"cell_type":"code","source":"!pip install hfda\nimport hfda","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:38:36.773047Z","iopub.execute_input":"2022-11-10T21:38:36.773625Z","iopub.status.idle":"2022-11-10T21:38:49.208169Z","shell.execute_reply.started":"2022-11-10T21:38:36.773589Z","shell.execute_reply":"2022-11-10T21:38:49.206805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"low pass filtering: <br>\n$ x_n = \\sum_{i=-k}^{+k} x_{n+i} $","metadata":{}},{"cell_type":"code","source":"def higuchi(x,k_max,winl=2048,overlap=0.5):\n    Nsamples = len(x)\n    nstart = int(winl*overlap)\n    nwindows = ((Nsamples - nstart)//winl)\n    win = np.arange(nstart, winl*nwindows, int(winl*overlap))\n    c=0\n    HFD = np.zeros(len(win))\n    for i in win:\n        signal = x[i-(winl//2):i+(winl//2)]\n        HFD[c] = hfda.measure(signal,k_max)\n        c+=1\n    \n    return win, HFD","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:38:49.210037Z","iopub.execute_input":"2022-11-10T21:38:49.210478Z","iopub.status.idle":"2022-11-10T21:38:49.219343Z","shell.execute_reply.started":"2022-11-10T21:38:49.210430Z","shell.execute_reply":"2022-11-10T21:38:49.218128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# By experiment we assume filtering window of length 3\nwin_sz = 3\nk_max = 10\n\nxticks = np.arange(0,600, 50)\n\nfor i in tqdm(range(1)):\n    z=np.random.randint(0,len(PREICTAL_PATHS)-1)\n    pre_mat = scipy.io.loadmat(PREICTAL_PATHS[z])\n    segment = dict(filter(lambda item: search_key in item[0], pre_mat.items()))\n    preictal_eeg = list(segment.values())[0][0][0][0]\n    mean_eeg = np.mean(preictal_eeg, 0)\n    # mean_eeg = np.convolve(mean_eeg, np.ones(win_sz)/win_sz, mode='same')\n    \n    win,h = higuchi(mean_eeg,18, winl=500)\n    freq,time,Z = scipy.signal.stft(mean_eeg,fs, nperseg=1000)\n    \n    fig0 = plt.figure(figsize=(30,5))\n    fig0.suptitle('MEAN DATA')\n    plt.plot(t0, mean_eeg)\n    plt.xticks(xticks)\n    plt.xlim((0,t0[-1]))\n    \n    fig1 = plt.figure(figsize=(30,5))\n    fig1.suptitle('HIGUCHI DIMENSION')\n    plt.xticks(xticks)\n    plt.plot(t0[win],h)\n    plt.xlim((0,t0[-1]))\n    plt.ylim((1,2.5))\n    \n    fig2 = plt.figure(figsize=(30,5))\n    fig2.suptitle('TIME-FREQUENCY')\n    plt.xticks(xticks)\n    plt.pcolormesh(time,freq,np.abs(Z)**2, vmin=0)\n    plt.ylim((0,50))\n    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-10T21:56:14.921322Z","iopub.execute_input":"2022-11-10T21:56:14.922367Z","iopub.status.idle":"2022-11-10T21:56:19.650464Z","shell.execute_reply.started":"2022-11-10T21:56:14.922326Z","shell.execute_reply":"2022-11-10T21:56:19.649345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in tqdm(range(1)):\n    z=np.random.randint(0,len(INTERICTAL_PATHS)-1)\n    inter_mat = scipy.io.loadmat(INTERICTAL_PATHS[z])\n    segment = dict(filter(lambda item: search_key in item[0], inter_mat.items()))\n    interictal_eeg = list(segment.values())[0][0][0][0]\n    mean_eeg = np.mean(interictal_eeg, 0)\n    # mean_eeg = np.convolve(mean_eeg, np.ones(win_sz)/win_sz, mode='same')\n    \n    win,h = higuchi(mean_eeg,18, winl=500)\n    freq,time,Z = scipy.signal.stft(mean_eeg,fs, nperseg=1000)\n    \n    fig0 = plt.figure(figsize=(30,5))\n    fig0.suptitle('MEAN DATA')\n    plt.plot(t0, mean_eeg)\n    plt.xticks(xticks)\n    plt.xlim((0,t0[-1]))\n    \n    fig1 = plt.figure(figsize=(30,5))\n    fig1.suptitle('HIGUCHI DIMENSION')\n    plt.xticks(xticks)\n    plt.plot(t0[win],h)\n    plt.xlim((0,t0[-1]))\n    plt.ylim((1,2.5))\n    \n    fig2 = plt.figure(figsize=(30,5))\n    fig2.suptitle('TIME-FREQUENCY')\n    plt.xticks(xticks)\n    plt.pcolormesh(time,freq,np.abs(Z)**2, vmin=0)\n    plt.ylim((0,50))\n    \nplt.show()","metadata":{},"execution_count":null,"outputs":[]}]}