<jupyter_start><jupyter_text>This notebook contains code to run the P300 Multi-subject EEG Guess The Number Dataset classification using a spiking neuralnetwork.To run the notebook, please download the dataset here: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/G9RRLNand place it in the **dataset** folder (do not rename it unless you change the name of the file in the notebook as well)The model used here is a slightly modified version of a CNN from the paper: https://www.sciencedirect.com/science/article/pii/S1746809419304185<jupyter_code># Necessary imports

import os
import nengo
import keras
import pandas as pd
from sklearn import metrics
import tensorflow as tf
import nengo_dl
import numpy as np
from scipy.io import loadmat
from sklearn.model_selection import train_test_split, ShuffleSplit
from tensorflow.python.keras import Input, Model
from tensorflow.python.keras.callbacks import EarlyStopping
from tensorflow.python.keras.layers import Conv2D, BatchNormalization, Dropout, AveragePooling2D, Flatten, Dense
from keras import backend as K
# First load the dataset and set seed for consistency
# The dataset should be by default saved in ../datasets/VarekaGTNEpochs.mat
# The file can be downloaded here https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/G9RRLN
dataset_path = os.path.join('dataset', 'VarekaGTNEpochs.mat')

params_folder_path = 'nengo_network_params' # path to saved parameters
os.makedirs('nengo_network_params', exist_ok=True)

seed = 0 # constant seed
np.random.seed(seed)
tf.random.set_seed(seed)
mat = loadmat(dataset_path)
target_data, non_target_data = mat['allTargetData'], mat['allNonTargetData'] # get target and non-target data
features = np.concatenate((target_data, non_target_data))

# Target labels are represented as (1, 0) vector, non target labels are represented as (0, 1) vector
target_labels = np.tile(np.array([1, 0]), (target_data.shape[0], 1)) # set 'target' as (1, 0) vector
non_target_labels = np.tile(np.array([0, 1]), (non_target_data.shape[0], 1)) # set 'non target' as (0, 1) vector
labels = np.vstack((target_labels, non_target_labels)) # concatenate target and non target labels

# Filter noise above 100 uV
threshold = 100.0
x_result, y_result = [], []
for i in range(features.shape[0]):
    if not np.max(np.abs(features[i])) > threshold:
        x_result.append(features[i])
        y_result.append(labels[i])

# Save data to numpy array
features, labels = np.array(x_result), np.array(y_result)
features = features.reshape((features.shape[0], 1, -1))
labels = labels.reshape((labels.shape[0], 1, -1))

# Print shapes
f'Features shape: {features.shape}, Labels shape: {labels.shape}'
# Define a function to create the CNN model
# This is a slightly modified version from https://www.sciencedirect.com/science/article/pii/S1746809419304185
def create_model():
    """
    Function to create tensorflow model
    """
    inp = Input(shape=(3, 1200, 1), name='input_layer')
    conv2d = Conv2D(filters=6, kernel_size=(3, 3), activation=tf.nn.relu)(inp)
    dropout1 = Dropout(0.5, seed=seed)(conv2d)
    avg_pooling = AveragePooling2D(pool_size=(1, 8), padding='same')(dropout1)
    flatten = Flatten()(avg_pooling)
    dense1 = Dense(100, activation=tf.nn.relu)(flatten)
    batch_norm = BatchNormalization()(dense1)
    dropout2 = Dropout(0.5, seed=seed)(batch_norm)
    output = Dense(2, activation=tf.nn.softmax, name='output_layer')(dropout2)

    return Model(inputs=inp, outputs=output)
# Since we used [1,0] as true (target) we also need to reference the correct index during metrics evaluation.
# In sklearn this means that we set the pos_label to 0 instead of 1
true_ref_idx = np.argmax(np.array([1, 0]))

def get_metrics(simulator, output_layer, x_test, y_test, minibatch_size, network_name):
    """
    Function for calculating metrics
    :param simulator: simulator instance
    :param input_layer: input layer reference
    :param output_layer: output layer reference
    :param x_test: features of the testing subset
    :param y_test: labels of the testing subset
    :param network_name: name of the network
    :return: accuracy, recall and precision metrics
    """

    # Truncate the remaining number of samples since the predict function uses minibatch_size as well
    samples = (x_test.shape[0] // minibatch_size ) * minibatch_size
    x_test, y_test = x_test[:samples], y_test[:samples]

    predictions = simulator.predict(x_test)[output_layer] # get result from the output layer when predicting on x_test
    predictions = predictions[:,-1,:] # get the last time step
    predictions_argm = np.argmax(predictions, axis=-1) # get the predicted labels

    y_test = np.squeeze(y_test, axis=1) # remove time dimension
    y_test_argm = np.argmax(y_test, axis=-1) # get labels, due to one-hot encoding 0 = target, 1 = non-target

    precision = metrics.precision_score(y_true=y_test_argm, y_pred=predictions_argm, pos_label=true_ref_idx) # get precision score
    recall = metrics.recall_score(y_true=y_test_argm, y_pred=predictions_argm, pos_label=true_ref_idx) # get recall
    f1 = metrics.f1_score(y_true=y_test_argm, y_pred=predictions_argm, pos_label=true_ref_idx) # get F1 score
    accuracy = metrics.accuracy_score(y_true=y_test_argm, y_pred=predictions_argm) # get accuracy

    # First column - "0" = target samples, second column "1" = non-target samples
    confusion_matrix = metrics.confusion_matrix(y_true=y_test_argm, y_pred=predictions_argm)

    # Log the statistics
    print(f'{network_name}: accuracy = {accuracy * 100}%, precision = {precision}, '
          f'recall = {recall}, f1 = {f1}')
    print('Confusion matrix:')
    print(confusion_matrix)

    return accuracy, precision, recall, f1, confusion_matrix

# Define functions to run the analog and spiking networks

def run_ann(model, train, valid, test, params_save_path, iteration, shuffle_training=True):
    """
    Run ANN via Nengo simulator. This fits the given model with the training data (train) and validates it using validation
    data (valid). Then metrics are calculated using the test data (test) and weights are saved to params_save_path
    :param shuffle_training: whether to shuffle data (default true)
    :param model: tensorflow model created from create_model() function
    :param train: pair of features and labels from training data
    :param valid: pair of features and labels from validation data
    :param test: pair of features and labels from test data
    :param params_save_path: output path to save weights of the network for SNN testing
    :return accuracy, precision, recall, f1 and confusion matrix from the testing data
    """

    # get features and labels for training, validation and testing data
    x_train, y_train = train[0], train[1]
    x_valid, y_valid = valid[0], valid[1]
    x_test, y_test = test[0], test[1]

    converter = nengo_dl.Converter(model)
    with nengo_dl.Simulator(converter.net, minibatch_size=64) as simulator:
        # Compile the model with binary cross-entropy and Adam optimizer
        simulator.compile(
            optimizer=keras.optimizers.Adam(),
            loss=keras.losses.BinaryCrossentropy(),
            metrics=['accuracy']
        )

        input_layer = converter.inputs[model.get_layer('input_layer')] # get nengo input layer
        output_layer = converter.outputs[model.get_layer('output_layer')] # get nengo output layer

        simulator.fit(
            x={ input_layer: x_train }, y={ output_layer: y_train },
            validation_data=({ input_layer: x_valid }, { output_layer: y_valid }),
            epochs=30,
            shuffle=shuffle_training,
            callbacks=[EarlyStopping(patience=5, verbose=1, restore_best_weights=True)] # early stop to avoid overfitting
        ) # train the model

        simulator.save_params(params_save_path) # save the params for the SNN

        accuracy, precision, recall, f1, confusion_matrix = get_metrics(simulator, output_layer, x_test, y_test,
                                                  minibatch_size=simulator.minibatch_size,
                                                  network_name=f'{iteration}. ANN')

        return {
            'accuracy': accuracy,
            'precision': precision,
            'recall': recall,
            'f1': f1,
            'confusion_matrix': confusion_matrix
        }


def run_snn(model, test, params_load_path, timesteps, scale_firing_rates, synapse, iteration):
    """
    Runs SNN on test data. Loads pre-trained weights from params_load path and uses timesteps, scale_firing_rates and synapse
    parameters for the spiking simulator.
    :param model: reference to the tensorflow model
    :param test: reference to the test features and labels
    :param params_load_path: path to the saved weights of the ANN
    :param timesteps: number of timesteps - i.e. how long is the input streamed to the network
    :param scale_firing_rates: firing rate scaling - increases the number of spikes
    :param synapse: synaptic smoothing
    :param iteration: iteration number (this is printed during the run)
    :return: accuracy, precision, recall, f1 and confusion matrix from the testing data
    """

    # Conversion of the TensorFlow model to a spiking Nengo model
    converter = nengo_dl.Converter(
        model=model,
        swap_activations={ tf.nn.relu: nengo.SpikingRectifiedLinear() },
        scale_firing_rates=scale_firing_rates,
        synapse=synapse
    )

    x_test, y_test = test[0], test[1] # get test features and labels

    with converter.net:
        nengo_dl.configure_settings(stateful=False)

    output_layer = converter.outputs[model.get_layer('output_layer')] # output layer for simulator
    x_test_time_tiled = np.tile(x_test, (1, timesteps, 1)) # tile x_test to match desired timesteps for simulator

    with nengo_dl.Simulator(converter.net, minibatch_size=41, progress_bar=False) as simulator:
        simulator.load_params(params_load_path)

        # Name of the network for print in get_metrics function
        name = f'{iteration}. SNN [timesteps={timesteps}, scale_firing_rates={scale_firing_rates}, synapse={synapse}]'
        accuracy, precision, recall, f1, confusion_matrix = get_metrics(simulator, output_layer, x_test_time_tiled, y_test,
                                                  minibatch_size=simulator.minibatch_size,
                                                  network_name=name)

        return {
            'accuracy': accuracy,
            'precision': precision,
            'recall': recall,
            'f1': f1,
            'confusion_matrix': confusion_matrix
        }
# Split the dataset into training and testing data
x_train, x_test, y_train, y_test = train_test_split(features, labels, test_size=0.25, random_state=seed, shuffle=True)

print(f'Train features shape: {x_train.shape}, train labels shape: {y_train.shape}')
print(f'Test features shape: {x_test.shape}, test labels shape: {y_test.shape}')

# Structures to store results from the ANN and the SNN
ann, snn = [], {}

# Configuration for the spiking network with format: timesteps, scale_firing_rates, synapse
snn_config = [
    [50, 1000, 0.01], # best performing parameters for simulator
    [50, 1000, None], # synaptic smoothing turned off
    [50, 1, 0.01], # spike scaling turned off
    [50, 1, None] # everything turned off, only RELU is swapped for spiking RELU
]

variants = []
for variant in snn_config:
        # name of the configuration
        name = f'snn [timesteps={variant[0]}, scaling={variant[1]}, synapse=None]' if variant[2] is None \
            else 'snn [timesteps={}, scaling={}, synapse={:3f}]'.format(variant[0], int(variant[1]), variant[2])
        variants.append(name)
        snn[name] = []

num_iterations = 30 # 30 iterations of CV
iteration = 1 # number of the current iteration
val_size = 0.25 # 25% of the data is used as validation data
for train_idx, val_idx in ShuffleSplit(n_splits=num_iterations, test_size=val_size, random_state=seed).split(x_train):
    print(f'Iteration: {iteration}')

    # Split all training data into current training data and validation data
    x_train_curr, y_train_curr = x_train[train_idx], y_train[train_idx]
    x_val_curr, y_val_curr = x_train[val_idx], y_train[val_idx]

    # Set params path
    params_path = os.path.join(params_folder_path, f'params_iter_{iteration}')

    # Create an untrained model
    model = create_model()

    # Run the analog network - train and evaluate
    ann_result = run_ann(model=model,
                         train=(x_train_curr, y_train_curr),
                         valid=(x_val_curr, y_val_curr),
                         test=(x_test, y_test),
                         params_save_path=params_path,
                         iteration=iteration
                         )
    K.clear_session() # clear session (otherwise PyCharm was sometimes causing memory leaks)
    ann.append(ann_result) # append the result to the ANN array

    # Create an array where results for each configuration of the spiking network will be stored
    snn_config_results = []
    for i, variant_conf in enumerate(snn_config):
        snn_result = run_snn(model=model,
                          test=(x_test, y_test),
                          params_load_path=params_path,
                          timesteps=variant_conf[0],
                          scale_firing_rates=variant_conf[1],
                          synapse=variant_conf[2],
                          iteration=iteration
                          )
        K.clear_session() # clear session
        snn[variants[i]].append(snn_result) # save results

    del model # delete the model (if this is not called it "may" create a memory leak - depends on the simulation machine)
    iteration += 1
# Create dictionary with the data for pandas dataframe
data = {
    'iterations': [x for x in range(1, num_iterations + 1)],
    'ann_accuracy': [x['accuracy'] for x in ann],
    'ann_precision': [x['precision'] for x in ann],
    'ann_recall': [x['recall'] for x in ann],
    'ann_f1': [x['f1'] for x in ann]
}

for variant_name in variants:
    data[f'{variant_name}_accuracy'] = [x['accuracy'] for x in snn[variant_name]]
    data[f'{variant_name}_precision'] = [x['precision'] for x in snn[variant_name]]
    data[f'{variant_name}_recall'] = [x['recall'] for x in snn[variant_name]]
    data[f'{variant_name}_f1'] = [x['f1'] for x in snn[variant_name]]
# Create pandas dataframe and save it to xlsx file
output_folder_path = 'p300_exp_output'
os.makedirs(output_folder_path, exist_ok=True)

df = pd.DataFrame(data)
df
df.to_excel(os.path.join(output_folder_path, 'data_iterations.xlsx'))

'File with iteration data successfully saved.'
network_names = ['ann'] + variants # names of each network in the pandas data frame
data_stats = {
    'models': network_names,
    'average_acc': [],
    'max_acc': [],
    'std_acc': [],
    'average_precision': [],
    'max_precision': [],
    'average_recall': [],
    'max_recall': [],
    'average_f1': [],
    'max_f1': []
} # statistics from the experiment

average_acc, max_acc, std_acc = [], [], []
average_precision, max_precision = [], []
average_recall, max_recall = [], []
for network_name in network_names:
    data_stats['average_acc'].append(df[f'{network_name}_accuracy'].mean())
    data_stats['max_acc'].append(df[f'{network_name}_accuracy'].max())
    data_stats['std_acc'].append(df[f'{network_name}_accuracy'].std())

    data_stats['average_precision'].append(df[f'{network_name}_precision'].mean())
    data_stats['max_precision'].append(df[f'{network_name}_precision'].max())

    data_stats['average_recall'].append(df[f'{network_name}_recall'].mean())
    data_stats['max_recall'].append(df[f'{network_name}_recall'].max())

    data_stats['average_f1'].append(df[f'{network_name}_f1'].mean())
    data_stats['max_f1'].append(df[f'{network_name}_f1'].max())
df_stats = pd.DataFrame(data_stats)

df_stats
df_stats.to_excel(os.path.join(output_folder_path, 'statistics.xlsx'))

'File with statistics successfully saved.'<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10263875123885034, 'fraction_numerical': 0.013875123885034688, 'mean_word_length': 3.7235225277940316, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '382388', 'license': 'MIT', 'n_tokens_mistral': 5166, 'n_tokens_neox': 4915, 'n_words': 1405, 'path': 'experiments/p300_experiment/p300_dataset_exp_convnet.ipynb', 'repo_name': 'honzikv/use-of-snn', 'size': '16144'}