<jupyter_start><jupyter_text>Train an LSTM classifier on fMRI activity of visual regions of the brain to predict type of visual stimuli Here, we will train a simple LSTM to predict which type of image was seen based on the brain activity of visual regions. For the first pass, we will simply classify whether the image came from ImageNet, COCO, or SUN. These databses have distinct types of images: focused single object, naturalistic objects, environmental scenes, respectively.We will train a separate LSTM for each of the ten brain regions that were extracted from each subject viewing these images. Perhaps we will find that in isolation neither region provides enough information to successfully classify, or one or two region will have particular higher accuracies. Whatever result we find will dictate how we proceed next. Maybe we can train combining all regions into one large training set, or use stacked LSTMs or hierarchical LSTMs. Let's find out using Keras.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import os
import pandas as pd

from keras.models import Sequential
from keras.layers import Dense, Dropout, LSTM, GRU
from sklearn.model_selection import train_test_split


data_dir = 'dataset/traindata/'<jupyter_output>Using TensorFlow backend.<jupyter_text>First pass Let's train an LSTM using some default hyperparameters.<jupyter_code>import pickle
from sklearn.feature_selection import SelectKBest, f_classif

"""
Load class numbers associated with each image as one-hot vectors
"""
def load_classes(Y, encoding = 'onehot'):
    if encoding not in ['onehot', 'interval']:
        raise ValueError('incorrect type of encoding')
    df = pd.read_csv(os.path.join(data_dir, 'image2labels.tsv'), sep='\t', header = None, names = ['ImgName', 'ClassNum', 'ClassLabel'])
    image2class = pd.Series(df.ClassNum.values,index=df.ImgName).to_dict()
    max_class_num = np.max(list(image2class.values()))
    if encoding == 'onehot':
        Y_out = np.zeros((Y.shape[0], max_class_num+2))
        for ind,img in enumerate(Y):
            if img == 'none':
                Y_out[ind][max_class_num+1] = 1
            else:
                Y_out[ind][image2class[img]] = 1
    elif encoding == 'interval':
        Y_out = np.zeros((Y.shape[0],))
        for ind,img in enumerate(Y):
            if img == 'none':
                Y_out[ind] = max_class_num + 1
            else:
                Y_out[ind] = image2class[img]
    
    return Y_out
        

"""
Subtract the mean signal from the resting state data
"""
def rest_normalize(X,Y):
    with open('dataset/traindata/sample2session.p', 'rb') as f:
        sample2session = pickle.load(f)
    subjects = list(set(sample2session['Subject']))
    baseline = {}
    for sub in subjects:
        rest = np.array([x for ind,x in enumerate(X) if Y[ind][16] == 1 and sample2session['Subject'].values[ind] == sub])
        baseline[sub] = np.mean(rest, axis = (0,1))
    X_new = np.array([x - baseline[sample2session['Subject'].values[ind]] for ind,x in enumerate(X)])
    
    return X_new    
    

"""
Load in the dataset.

Return train, dev, test splits, standardized
"""
def load_data(dim_reduce = False, num_dims = 1000):
    # Load from disk
    with open(data_dir + 'X_raw.p', 'rb') as f:
        X = np.array(pickle.load(f))
    with open(data_dir + 'Y_raw.p', 'rb') as f:
        Yimages = np.array(pickle.load(f))
    if dim_reduce:
        classes = load_classes(Yimages, encoding = 'interval')
        X_mean = np.mean(X, axis=1)
        kbest = SelectKBest(f_classif, k = num_dims).fit(X_mean,classes)
        best_voxels = np.argsort(kbest.scores_)[::-1][:num_dims]
        X = X[:,:,best_voxels]
    # Convert class labels as one hot vectors
    classes = load_classes(Yimages)
    # Standardize data
    X_scaled = rest_normalize(X, classes)
    # Split into train,dev,test
    X_train, X_devtest, Y_train, Y_devtest = train_test_split(X_scaled, classes, test_size = 0.4)
    X_dev, X_test, Y_dev, Y_test = train_test_split(X_devtest, Y_devtest, test_size = 0.5)
    
    print('%20s %20s' % ('Shape of X_train:', str(X_train.shape)))
    print('%20s %20s' % ('Shape of Y_train:', str(Y_train.shape)))
    print('%20s %20s' % ('Shape of X_dev:', str(X_dev.shape)))
    print('%20s %20s' % ('Shape of Y_dev:', str(Y_dev.shape)))
    print('%20s %20s' % ('Shape of X_test:', str(X_test.shape)))
    print('%20s %20s' % ('Shape of Y_test:', str(Y_test.shape)))
    
    return X_train, X_dev, X_test, Y_train, Y_dev, Y_test
# Load with reduced dimensions
X_train, X_dev, X_test, Y_train, Y_dev, Y_test = load_data(dim_reduce = False)
"""
Create LSTM model using given hyperparameters

"""
def create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = 0.5):
    model = Sequential()
    model.add(GRU(n_hidden, input_shape=(n_timesteps,n_features)))
    model.add(Dropout(dropout))
    #model.add(Dense(n_hidden, activation='relu'))
    model.add(Dense(n_classes, activation='softmax'))
    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    
    return model
n_timesteps, n_features, n_classes = X_train.shape[1], X_train.shape[2], Y_train.shape[1]
n_hidden = 100
dropout = 0.5

model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
model.summary()
"""
Train model and output train and dev accuracy history

"""

def train_model(model, X_train, X_dev, Y_train, Y_dev, epochs, batch_size = 64, verbose = 1):
    train_history = model.fit(X_train, Y_train, validation_data = (X_dev, Y_dev), epochs=epochs, batch_size=batch_size, verbose=verbose)
    train_acc,dev_acc = train_history.history['accuracy'], train_history.history['val_accuracy']
    
    return model,train_acc, dev_acc
model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, 30, batch_size = 64)
"""
Plot the training history
"""
def plot_acc(train_acc, dev_acc, title = 'LSTM model accuracy'):
    plt.close('all')
    fig, ax = plt.subplots(1,1, figsize = (6,4))
    ax.plot(train_acc, color = 'blue', label = 'train')
    ax.plot(dev_acc, color = 'green', label = 'dev')
    ax.set_xlabel('Epoch')
    ax.set_ylabel('Accuracy')
    ax.set_title(title)
    ax.legend()
    plt.show()
plot_acc(train_acc, dev_acc, title = 'GRU model accuracy')
dropout = 0.5
n_hidden = 500
model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, epochs = 20, batch_size = 64, verbose = 1)<jupyter_output>Train on 11628 samples, validate on 3876 samples
Epoch 1/20
 4480/11628 [==========>...................] - ETA: 1:46 - loss: 3.5481 - accuracy: 0.0879<jupyter_text>It is clearly overfitting, maybe we can try to reduce the variance. How about increasing dropout?<jupyter_code>dropout = 0.8
n_hidden = 100
model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, epochs = 20, batch_size = 64, verbose = 1)
plot_acc(train_acc, dev_acc)<jupyter_output><empty_output><jupyter_text>That did not improve our validation error. How about keeping dropout but changing the number of hidden nodes?<jupyter_code>dropout = 0.5
n_hidden = 50
model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, epochs = 30, batch_size = 64, verbose = 1)
plot_acc(train_acc, dev_acc)<jupyter_output><empty_output><jupyter_text>Decreasing number of hidden units may not help either... the other alternative is to get more data. We may have to use multiple regions to train one model. We can first train an individual model for each region to see if any region can get good accuracy. Then we can try combining the regions. If none of those improve performance, then we'll have to look into combining multiple LSTMs.  One model per region<jupyter_code>n_hidden = 100
dropout = 0.5
epochs = 50

masks = set([f.split('_')[1][:-4] for f in os.listdir(data_dir) if os.path.isfile(os.path.join(data_dir,f)) and '.npy' in f and 'all' not in f])

region_acc = {}

for mask in masks:   
    X_train, X_dev, X_test, Y_train, Y_dev, Y_test = load_region_data(mask)
    n_timesteps, n_features, n_classes = X_train.shape[1], X_train.shape[2], Y_train.shape[1]

    model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
    model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, epochs = epochs, batch_size = 64, verbose = 0)
    
    region_acc[mask] = (train_acc[-1], dev_acc[-1])
    
    plot_acc(train_acc, dev_acc, title = mask)<jupyter_output>Shape of X_train:      (11628, 5, 100)
   Shape of Y_train:          (11628, 17)
     Shape of X_dev:       (3876, 5, 100)
     Shape of Y_dev:           (3876, 17)
    Shape of X_test:       (3876, 5, 100)
    Shape of Y_test:           (3876, 17)<jupyter_text>No region looks promising... maybe we could try combining the regional data. Just worried about curse of dimensionality.  One model for all regions<jupyter_code>"""
Load in all the regions
"""
def load_all_data():
    
    masks = set([f.split('_')[1][:-4] for f in os.listdir(data_dir) if os.path.isfile(os.path.join(data_dir,f)) and '.npy' in f and 'all' not in f])

    Y_all = np.load(os.path.join(data_dir,'Ylabels_LHLOC.npy'))
    n_samples = Y_all.shape[0]

    X_all = [np.load(os.path.join(data_dir,'X_' + region + '.npy'))[:n_samples] for region in masks]
    X_all = np.concatenate(X_all, axis = 2)
    
    classes = load_classes(Y_all)
    
    # Standardize data
    X_scaled = np.array([(x - np.mean(x,axis=1)[:,np.newaxis]) / np.std(x,axis=1)[:,np.newaxis] for x in X_all])
    # Split into train,dev,test
    X_train, X_devtest, Y_train, Y_devtest = train_test_split(X_scaled, classes, test_size = 0.4)
    X_dev, X_test, Y_dev, Y_test = train_test_split(X_devtest, Y_devtest, test_size = 0.5)
    
    print('%20s %20s' % ('Shape of X_train:', str(X_train.shape)))
    print('%20s %20s' % ('Shape of Y_train:', str(Y_train.shape)))
    print('%20s %20s' % ('Shape of X_dev:', str(X_dev.shape)))
    print('%20s %20s' % ('Shape of Y_dev:', str(Y_dev.shape)))
    print('%20s %20s' % ('Shape of X_test:', str(X_test.shape)))
    print('%20s %20s' % ('Shape of Y_test:', str(Y_test.shape)))
    
    return X_train, X_dev, X_test, Y_train, Y_dev, Y_test
n_hidden = 100
dropout = 0.5
epochs = 50

X_train, X_dev, X_test, Y_train, Y_dev, Y_test = load_all_data()

n_timesteps, n_features, n_classes = X_train.shape[1], X_train.shape[2], Y_train.shape[1]

model = create_model(n_timesteps, n_features, n_hidden, n_classes, dropout = dropout)
model, train_acc, dev_acc = train_model(model, X_train, X_dev, Y_train, Y_dev, epochs = epochs, batch_size = 64, verbose = 1)
plot_acc(train_acc, dev_acc)<jupyter_output><empty_output><jupyter_text>Yikes, no improvement... perhaps we should reconsider our architecture or try to find more data. What classes are we getting wrong?<jupyter_code>from sklearn.metrics import confusion_matrix

Y_pred = model.predict(X_test)
matrix = confusion_matrix(Y_test.argmax(axis=1), Y_pred.argmax(axis=1))
print(matrix)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12411410347271438, 'fraction_numerical': 0.022147413182140325, 'mean_word_length': 4.197513812154696, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '300598', 'license': 'MIT', 'n_tokens_mistral': 4156, 'n_tokens_neox': 3951, 'n_words': 1100, 'path': 'LSTM.ipynb', 'repo_name': 'arashjamalian/fmriNet', 'size': '11288'}