# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Use machine learning modeling to investigate which gamma frequency band represents the contrast level of grating stimuli.
#
# This notebook demonstrates the pipeline of using __Support Vector Machine (SVM)__ to test which gamma frequency band can predict the contrast information of grating stimulus presented to the patient subject.
#
# Check the following link for the research story of why I am doing it: https://www.cell.com/current-biology/fulltext/S0960-9822(19)31020-6?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0960982219310206%3Fshowall%3Dtrue#secsectitle0035
#
# The following analysis is relevant to the result session __"Gratings Contrast Can Be Decoded Better Using NBG Than BBG"__, and method session __"Visual grating contrast classification"__.
#
# Two frequency band to compare are: 20-60 Hz (variable name in the code is 'NBG'), 70-150 Hz (variable name in the code is 'BBG').
#
# This notebook will: 
# <br>
# 1) read event files to extract contrast information of trials from the task.
# <br>
# 2) read feature file to extract feature matrix as the input for the SVM modeling.
#
# Note: 
# <br>
# 1) This notebook won't provide the codes for the preprocessing step for the ECoG data. The data is already clean and transform into time-frequency data using Morlet wavelet transform (codes were not provided in the repositories). You can take a loot at 'extract_epoch.m' (use in MATLAB) to see how I extract the task epoch into the feature file I need for the modeling. However, you may have a different way to structure your EEG/ECoG file and may do it in a different way.
# <br>
# 2) This notebook mainly demonstrate how the modeling is done in one single electrode in one patients.

import os
import scipy.io as sio
import numpy as np
import pandas as pd
import itertools
from sklearn import svm,preprocessing
from sklearn.model_selection import LeaveOneOut, permutation_test_score
from sklearn.metrics import accuracy_score
from sklearn.pipeline import Pipeline

# ### 1. Read event file to extract the condition labels of task trials

# Set the path of the file and read the mat file of condition events. I have provided one example event file from one task block of one patient subject in this repository.

# +
# set path and read event file
# modify the event_path based on how you save your data

event_path   = 'example_data/task_events_example.mat'
event_file   = sio.loadmat(event_path)
# -

# In this project, there are trials of a given condition that I want to exclude, here I delete the events of those trials (trials present patients with a grating stimulus at the orientation of 45 degree).

# get contrast labels and detele those at 45 degree
contrast     = np.asscalar(event_file['events_info']['trial_contrast'])
odd_orient   = np.asscalar(event_file['events_info']['trial_orientation']) 
contrast_new = np.delete(contrast, np.where(odd_orient == 45)[0]) 

#  I recode the task condition into 1,2,3 for three different task conditions.

# label the task conditions (1 = 20% contrast, 2 = 50% contrast, 3 = 100% contrast)
label_v                 = np.reshape(contrast_new,(90,))
label_v[label_v == 1]   = 3
label_v[label_v == 0.5] = 2
label_v[label_v == 0.2] = 1

# Cool, the label_v is the label file we want for the machine learning training process (the 'observed values' of output in a prediction question).
#
# Next, I will extract the feature matrix for the model training and prediction.

# ### 2. Read the time-series data
# With this notebook, the real data of one electrode from our patients is provided in the repository.

# read time series data
# the data with be the spectrogram data of one task block
tf_path = 'example_data/tf_epochs_example.mat'
tf_data = sio.loadmat(tf_path) # it is a dict file, the "tf_epochs" in it is the features we need

# Take a look at the shape of the data to understand it.

tf_data['tf_epochs'].shape

# Dimensions of tf_data['tf_epochs'] mean:
# <br>
# "200" frequency points, from 2 Hz to 201 Hz;
# <br>
# "250" time poins, from 250 ms to 500 ms after the stimulus onset;
# <br>
# "90" trials of the task -- grating stimuli at three contrast levels were shown to the subject for 90 times in one task.
#
# We need to take the average of the data across the time window:

tf_avg = np.mean(tf_data['tf_epochs'],axis =1)

# Next, define the variable of the boundary of two gamma frequency bands and extract the feature input matrixs.

# +
# define two gammas
NBG = [20,60]
BBG = [70,150]

# get feature inputs
feature_NBG = tf_avg[(NBG[0]-2):(NBG[1]-1),:].T # transform the matrix to make the column is frequency info
feature_BBG = tf_avg[(BBG[0]-2):(BBG[1]-1),:].T


# -

# Finally we can run the model.
#
# The model pipeline will:
# <br>
# 1) Scale the data using MinMax method.
# <br>
# 2) Implement a leave-one-out-cross-validation method to get the prediction accuracy.
# <br>
# 3) Get the significant level of the prediction accuracy using permutation test. (In this notebook we will use 1000-time-permutation).
#
# In the brain science, the aim is not to create a fancy model to have the highest prediction performance, but to use modeling method to compare the prediction performance of different brain signals. In this way, we can figure out whether a given brain function is represented by a given brain signal, or the brain signal from a given brain region.
#
# Here we want to compare:
# <br>
# 1) The prediction performance of model using NBG as input vs using BBG as input -- so we can see which frequency band is representing the stimulus contrast information.
# <br>
# 2) The prediction performance of model using NBG within visual cortex vs outside of visual cortex -- it serves as a great control analysis to validate the SVM model, the NBG from areas outside of visual cortex should not be able to predict the stimulus visual properties.
#

# function of the SVM modeling
def linear_SVM(feature_mx,npermu,label):
    """This function conduct SVM on single electrode data.
    
    feature_mx: the feature input of the model.
    
    nperm: number of permutation times. 
    
    label: label vector. """
 
    # SVM pipeline
    loo = LeaveOneOut()
    pipe = Pipeline([('scaler', preprocessing.MinMaxScaler()),('clf', svm.SVC(kernel='linear'))])
    score, permutation_scores, pvalue = permutation_test_score(pipe, feature_mx, label, scoring = 'accuracy', cv=loo, n_permutations=npermu, n_jobs = -1)

    return score,pvalue


# get the prediction results
score_NBG, pvalue_NBG = linear_SVM(feature_NBG,1000,label_v)
score_BBG, pvalue_BBG = linear_SVM(feature_BBG,1000,label_v)
print('the accuracy of model prediction using NBG is:' + str(score_NBG) + '; p-value: ' + str(pvalue_NBG))
print('the accuracy of model prediction using BBG is:' + str(score_BBG) + '; p-value: ' + str(pvalue_BBG))

# From this example electrode we can see the NBG is representing the stimulus visual property -- contrast level. 
#
# Next I will extract the data of an electrode from areas outside of visual cortex and do the same modeling again.
#

# read data and make input feature
tf_path = 'example_data/tf_epochs_example_nonvisual.mat'
tf_data = sio.loadmat(tf_path)
tf_avg = np.mean(tf_data['tf_epochs'],axis =1)
feature_nonvisual = tf_avg[(NBG[0]-2):(NBG[1]-1),:].T

# get the prediction result
score_nonvisual, pvalue_nonvisual = linear_SVM(feature_nonvisual,1000,label_v)
print('the accuracy of model prediction using NBG outside of visual area is:' + str(score_nonvisual) + '; p-value: ' + str(pvalue_nonvisual))

# As expected, the NBG outside of visual area cannot predict the contrast level of the visual stimulus well, so we know the model works.
# Next, loop through all the visual electrodes from all the patients to get the prediction accuracy and p-value for each electrode.

# ### 3. Loop through data from all electrodes to the SVM pipeline
# Due to the big size of our raw data, I can't provide all the data we used for this project in the repository. However, I will demonstrate the way I save the subject information in the dataframe and the code to loop through the data using the dataframe.
#
# Let me first make a "fake" subject information sheet first.

# create the datafram that saves the subject information
d = {'subject_code': ['AAA','AAB','AAC','AAD'],
     'task_block_code': ['001,013','009','014,023','005,006'],
     'electrode_code': ['097,098,099,100,122,123','065,066,077,088,089','001,002,003','123,124,125,126,127']}
df = pd.DataFrame(data = d)

df

# In my real project, I have more subjects, more task blocks and more electrodes. In addition, I separate the electrodes into "visual" and "non_visual", here I just made the dataframe simple enough to understand.

# create a dataframe to save the result
# I need to record: subject code, electrode_code, accuracy, p-value
# Note: for a given electrode, data from multiple task blocks will be concatenate together for SVM modeling
d_r = {'subject','elec','acc','p'}
df_result = pd.DataFrame(columns=d_r)

df_result



StarcoderdataJupyterScriptsDedupFiltered