<jupyter_start><jupyter_text>RSA assignment <NAME>  part 1 - loading data: First we import the data from mne.datasets and save it on our computer. This is ~6GB so don't run it. We will loosely follow: https://mne.tools/dev/auto_examples/decoding/decoding_rsa.htmlsphx-glr-auto-examples-decoding-decoding-rsa-py in this first part of the assignment.<jupyter_code>## don't run:
import mne
from mne.datasets import visual_92_categories
data_path = visual_92_categories.data_path()<jupyter_output>Using default location ~/mne_data for visual_92_categories...
Downloading archive MNE-visual_92_categories-data-part1.tar.gz to C:\Users\95\mne_data
Downloading https://files.osf.io/v1/resources/rxvq7/providers/osfstorage/59c600249ad5a1025db41c83?version=1&action=download&direct (3.47 GB)<jupyter_text>importing packages & tools Explanation for usage in comments.<jupyter_code>import os.path as op ## needed to paste variables together with strings (to get strings for data-paths). 
import numpy as np ## numpy needed in several places for working with arrays. 
from pandas import read_csv ## needed to read the data that we have stored. 
import matplotlib.pyplot as plt ## pyplot needed for plotting later. 
from mne.io import read_raw_fif, concatenate_raws ## reading certaing file-types.<jupyter_output><empty_output><jupyter_text>Reading data Reading the meta-data (stimuli-trigger mapping) data from csv & restricting the number of conditions for computational feasibility.<jupyter_code>## reading data. 
fname = op.join(data_path, 'visual_stimuli.csv') ## path & file concatenated to one string. 
conds = read_csv(fname) ## read the csv. 
print(conds.head(5)) ## checking the data
print(conds.tail(5)) ## checking the data. 

## restrict number of conditions: 
max_trigger = 24
conds = conds[:max_trigger]  # take only the first 24 rows
for c in conds.values: 
    print(c[0:2])<jupyter_output>[0 'human bodypart']
[1 'human bodypart']
[2 'human bodypart']
[3 'human bodypart']
[4 'human bodypart']
[5 'human bodypart']
[6 'human bodypart']
[7 'human bodypart']
[8 'human bodypart']
[9 'human bodypart']
[10 'human bodypart']
[11 'human bodypart']
[12 'human face']
[13 'human face']
[14 'human face']
[15 'human face']
[16 'human face']
[17 'human face']
[18 'human face']
[19 'human face']
[20 'human face']
[21 'human face']
[22 'human face']
[23 'human face']<jupyter_text>Converting the format We want a list with strings reflecting the format of each of the conditions.<jupyter_code>## stimulus trigger mapping: 
conditions = [] ## empty list. 
for c in conds.values: ## looping over the values in the pandas dataframe. 
    cond_tags = list(c[:2]) ## first two columns of each row (0-2)
    cond_tags += [('not-' if i == 0 else '') + conds.columns[k]
                  for k, i in enumerate(c[2:], 2)] ## adds "not" when there is 0 in column.
    conditions.append('/'.join(map(str, cond_tags)))
print(conditions[:10])
## turn them into a dictionary: 
event_id = dict(zip(conditions, conds.trigger + 1))
event_id['0/human bodypart/human/not-face/animal/natural']<jupyter_output><empty_output><jupyter_text>Reading the actual MEG data (fif)<jupyter_code>## read MEG data 
n_runs = 4 # 4 for full data-set. 
fname = op.join(data_path, 'sample_subject_%i_tsss_mc.fif') # same trick with op.join. 
raws = [read_raw_fif(fname % block, verbose='error')  # ignore filename warnings
        for block in range(n_runs)] 
raw = concatenate_raws(raws) 

events = mne.find_events(raw, min_duration=.002) 

events = events[events[:, 2] <= max_trigger]
## which picks do we want?
picks = mne.pick_types(raw.info, 
                       meg=True) ## MEG-data 

## picking out epochs. 
epochs = mne.Epochs(raw, 
                    events=events, 
                    event_id=event_id, ## here we have specified our IDs. 
                    baseline=None,
                    picks=picks, ## defined just before. 
                    tmin=-.1, ## 100 ms before trigger. 
                    tmax=.5, ## 500 ms after trigger. 
                    preload=True)

## looking at the data: 
epochs['face'].average().plot();
epochs['not-face'].average().plot();<jupyter_output><empty_output><jupyter_text>Downsampling<jupyter_code>## downsample epochs from 1000 hz to 200 hz. 
epochs.resample(sfreq = 200, n_jobs = -1) #200 hz.<jupyter_output>[Parallel(n_jobs=8)]: Using backend LokyBackend with 8 concurrent workers.
[Parallel(n_jobs=8)]: Done  10 tasks      | elapsed:    3.7s
[Parallel(n_jobs=8)]: Done 272 tasks      | elapsed:    3.8s
[Parallel(n_jobs=8)]: Done 33808 tasks      | elapsed:    6.4s
[Parallel(n_jobs=8)]: Done 98320 tasks      | elapsed:   11.1s
[Parallel(n_jobs=8)]: Done 181264 tasks      | elapsed:   17.3s
[Parallel(n_jobs=8)]: Done 220320 out of 220320 | elapsed:   19.9s finished<jupyter_text>Creating results matrixHas to have dimension 24 x 24 x number of time-points. This is because we have 24 unique pictures that we want to predict. This means that our results matrix will have all possible combinations of those 24 pictures (x and y) and we have this for each time-point.<jupyter_code>## creating results matrix (numpy array) of dimension 24x24xtime
results = np.zeros([24, 24, epochs.get_data().shape[-1]]) #or epochs.get_data().shape[-1]
results.shape<jupyter_output><empty_output><jupyter_text>Importing from sklearn Explanation in comments.<jupyter_code>from sklearn.model_selection import StratifiedKFold ## for cross-validation
from sklearn.pipeline import make_pipeline ## pipeline for 
from sklearn.preprocessing import StandardScaler ## scaling variables after train/test split. 
from sklearn.linear_model import LogisticRegression ## logistic regression is our ML-algorithm of choice. 
from sklearn.metrics import roc_auc_score ## I don't think we actually use this in the end..
from mne.decoding import SlidingEstimator, cross_val_multiscore ## sliding estimator wraps pipeline & cross-validation.<jupyter_output><empty_output><jupyter_text>preparing the pipeline & modelFirst we make a pipeline which first scales the data and then applies logistic regression.<jupyter_code>clf = make_pipeline(StandardScaler(), # computed as: z = (x - u) / sd [x = value, u = mean, sd = standard deviation]
                      LogisticRegression()) ## using baseline settings. (l2-norm and C = 1)

time_decoder = SlidingEstimator(
    clf, #the pipeline. 
    n_jobs = -1, ## using all our cores.
    scoring='roc_auc') ## using area under the curve as our performance metric. 

cv = StratifiedKFold() #should have randomstate (but library.config issues) - do we even use this?<jupyter_output><empty_output><jupyter_text>Running everything What we want:  We want to to train logistic regression on each unique combination of two pictures. To use all of our data we do this by cross-validation. This means that a lot of cross-validation runs will be run. We get a score (roc-auc) for each cross-validation fold (we have 5). We then take the mean of this score and put into the (initially empty) results This has to be done for each time-point (what the time_decoder takes care of) and - obviously - it has to be done for each (unique) combination of pictures. Thus, it takes a lot of time. I have already run this and saved the file as "results_4run.npy".<jupyter_code>for ii in range(24): ## 24 pictures
    for jj in range(24): ## 24 pictures
        if ii > jj: ## Only do something if this is the case. Ensures only below diag. 
            cond_0 = ii ## define cond. 0 (int between 0-23)
            cond_1 = jj ## define cond. 1 (int. between 0-23)
            X = epochs[str(cond_0), str(cond_1)].get_data() ## taking X out. 
            y = epochs[str(cond_0), str(cond_1)].events[:, 2] ## trigger values (each diff cat)
            score = cross_val_multiscore(time_decoder, ## across all time-points. 
                                         X, ## train on two conditions (pictures)
                                         y, ## test on the same two conditions (pictures)
                                         cv = 5, ## five fold cross-validation
                                        n_jobs = -1) ## important: use all cores. 
            results[cond_0, cond_1] = score.mean(axis = 0) ## mean across folds axis and input to results array.<jupyter_output>[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   5 | elapsed:    6.3s remaining:    9.4s
[Parallel(n_jobs=4)]: Done   5 out of   5 | elapsed:    9.3s finished
[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   5 | elapsed:    5.8s remaining:    8.7s
[Parallel(n_jobs=4)]: Done   5 out of   5 | elapsed:    9.0s finished
[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   5 | elapsed:    5.6s remaining:    8.4s
[Parallel(n_jobs=4)]: Done   5 out of   5 | elapsed:    8.6s finished
[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   5 | elapsed:    5.6s remaining:    8.4s
[Parallel(n_jobs=4)]: Done   5 out of   5 | elapsed:    8.6s finished
[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Don[...]<jupyter_text>Save the results for later<jupyter_code>## np.save("results_4run.npy", results) ## save stuff. 
results = np.load("results_4run.npy")<jupyter_output><empty_output><jupyter_text>Plotting 1First we would like to plot our average accuracy across time. I.e., we will average over all the comparisons for each time-point.  The output (run below) is not impressive. We only just manage to get accuracy > 70% around 100 ms. after stimulus on-set.<jupyter_code>## accuracy across time: 
tri_idx = np.tril_indices(24, k=-1) #all different combinations of lower triangle. k-1 (no diag.)
plt.plot(epochs.times, results[tri_idx[0], tri_idx[1]].mean(axis = 0)) # mean over comparisons axis.<jupyter_output><empty_output><jupyter_text>Plotting 2Now let us plot accuracy for some different time-points without averaging the accuracies at pair-wise classification. We will do: * -100 ms.* 100 ms. * 300 ms. * 500 ms.recall that sampling frequency is 200 hz. First we define time-points and time-labels and compute the corresponding images.<jupyter_code>timePoints = [0, 40, 80, 119]
timeLabels = ["-100 ms", "100 ms", "300 ms", "500 ms"]
images = [results[:, :, x] for x in timePoints]<jupyter_output><empty_output><jupyter_text>Then we define a function to nicely show these images. This was inspired by some online code.<jupyter_code>## function (half stolen from the internet): 
def show_images(images, cols = 1, titles = None):
    """Display a list of images in a single figure with matplotlib.
    
    Parameters
    ---------
    images: List of np.arrays compatible with plt.imshow.
    
    cols (Default = 1): Number of columns in figure (number of rows is 
                        set to np.ceil(n_images/float(cols))).
    
    titles: List of titles corresponding to each image. Must have
            the same length as titles.
    """
    
    assert((titles is None) or (len(images) == len(titles)))
    
    n_images = len(images)
    if titles is None: titles = ['Image (%d)' % i for i in range(1,n_images + 1)]
        
    fig = plt.figure()
    for n, (image, title) in enumerate(zip(images, titles)):
        a = fig.add_subplot(cols, np.ceil(n_images/float(cols)), n + 1)
        plt.imshow(image, cmap = 'viridis')
        plt.axline((11.5, 0), (11.5, 23), color = "black")
        plt.axline((0, 11.5), (23, 11.5), color = "black")
        plt.colorbar() 
        a.set_title(title)
    fig.set_size_inches(np.array(fig.get_size_inches()) * n_images)
    plt.show()<jupyter_output><empty_output><jupyter_text>Then we run the function on our images and labels. We also specify cols = 2 to get a nice 2-by-2 grid presentation. It is actually clear from the picture that accuracy is best around ~100-300 ms.<jupyter_code>## running the function to generate plots: 
show_images(images = images, cols = 2, titles = timeLabels)<jupyter_output><ipython-input-38-0174263a7145>:23: MatplotlibDeprecationWarning: Passing non-integers as three-element position specification is deprecated since 3.3 and will be removed two minor releases later.
  a = fig.add_subplot(cols, np.ceil(n_images/float(cols)), n + 1)<jupyter_text>plotting 3Plotting accuracy per group over time.groups selected: 1. within not-face (should be hard - low accuracy)2. within is-face (should be hard - low accuracy)3. across condition (should be easiest - high accuracy)<jupyter_code>epochs # 0-11 (not face), 12-23 (face)
not_face = results[0:12,0:12,:] #smaller triangular (upper left).
is_face = results[12:24, 12:24, :] #smaller triangular (lower right).
across = results[11:23, 0:12] #across conditions (lower left). 
tri_idx = np.tril_indices(12, k=-1) #can be used for (1) and (2)
square_idx = np.indices((12, 12)) #used for condition (3)
## plotting them: 
plt.plot(epochs.times, not_face[tri_idx[0], tri_idx[1]].mean(axis = 0),
        label = "not faces")
plt.plot(epochs.times, is_face[tri_idx[0], tri_idx[1]].mean(axis = 0),
        label = "faces")
plt.plot(epochs.times, across[square_idx[0].flatten(), square_idx[1].flatten()].mean(axis = 0),
        label = "across conditions")
plt.legend()
plt.title('AUC across conditions')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11309965839893064, 'fraction_numerical': 0.03267488489529185, 'mean_word_length': 3.787415570565233, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 64, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '365203', 'license': 'MIT', 'n_tokens_mistral': 4630, 'n_tokens_neox': 4253, 'n_words': 1450, 'path': 'practical_classes/week_8/RSA.ipynb', 'repo_name': 'victor-m-p/adv_cognitive_neuroscience', 'size': '13473'}