eeg_biometrics_autoencoder.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Autoencoder",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/BRomans/IdMind/blob/main/eeg_biometrics_autoencoder.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1aTCrOSwlbGd"
      },
      "source": [
        "## Load libraries & initialise environment"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8sQFvwJpfLSt"
      },
      "source": [
        "# import libraries\r\n",
        "import os\r\n",
        "import csv\r\n",
        "import random\r\n",
        "import numpy as np\r\n",
        "import pandas as pd\r\n",
        "import tensorflow as tf\r\n",
        "import matplotlib.pyplot as plt\r\n",
        "import seaborn as sns\r\n",
        "sns.set(color_codes=True)\r\n",
        "%matplotlib inline\r\n",
        "\r\n",
        "from google.colab import drive\r\n",
        "from sklearn.preprocessing import MinMaxScaler, RobustScaler, QuantileTransformer\r\n",
        "from copy import deepcopy\r\n",
        "from sklearn.externals import joblib\r\n",
        "from numpy.random import seed\r\n",
        "from keras.layers import Input, Dropout, Dense, LSTM, TimeDistributed, Reshape, \\\r\n",
        "          RepeatVector, MaxPooling1D, Conv1D, Flatten, Conv1DTranspose, UpSampling1D, \\\r\n",
        "          AveragePooling1D\r\n",
        "from keras.models import Model, Sequential\r\n",
        "from keras.optimizers import Adam\r\n",
        "from keras import regularizers\r\n",
        "from keras import backend as K\r\n",
        "from keras.callbacks import ModelCheckpoint, EarlyStopping\r\n",
        "from keras.utils import to_categorical\r\n",
        "\r\n",
        "from tensorflow.python.client import device_lib\r\n",
        "# print(device_lib.list_local_devices())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gXQT6VF9fRAK"
      },
      "source": [
        "drive.mount(\"/content/drive\")\r\n",
        "dirpath = \"/content/drive/MyDrive/ml2-eeg-biometrics/train-test-data/\" "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E6YB2XBaXlTy"
      },
      "source": [
        "# Seed value\r\n",
        "# Apparently you may use different seed values at each stage\r\n",
        "seed_value = 10\r\n",
        "tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR)\r\n",
        "\r\n",
        "# 1. Set `PYTHONHASHSEED` environment variable at a fixed value\r\n",
        "os.environ['PYTHONHASHSEED']=str(seed_value)\r\n",
        "\r\n",
        "# 2. Set `python` built-in pseudo-random generator at a fixed value\r\n",
        "\r\n",
        "random.seed(seed_value)\r\n",
        "\r\n",
        "# 3. Set `numpy` pseudo-random generator at a fixed value\r\n",
        "np.random.seed(seed_value)\r\n",
        "\r\n",
        "# 4. Set the `tensorflow` pseudo-random generator at a fixed value\r\n",
        "tf.compat.v1.set_random_seed(seed_value)\r\n",
        "\r\n",
        "# 5. Configure a new global `tensorflow` session\r\n",
        "session_conf = tf.compat.v1.ConfigProto(intra_op_parallelism_threads=1, inter_op_parallelism_threads=1)\r\n",
        "sess = tf.compat.v1.Session(graph=tf.compat.v1.get_default_graph(), config=session_conf)\r\n",
        "tf.compat.v1.keras.backend.set_session(sess)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SeTPxR0aYNxV"
      },
      "source": [
        "## Load & Process Data\r\n",
        "\r\n",
        "\r\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BInAGYqORGnb"
      },
      "source": [
        "##### Load data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vlN5H_lIhPin"
      },
      "source": [
        "To begin with, we load the 3 parts of the dataset, training, test and validation that we split in the pre-processing phase."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EAqlAto997Wp"
      },
      "source": [
        "def load_datasets():\r\n",
        "  x_train = np.asarray(np.load(dirpath + 'x_train.npy')).astype(np.float32)\r\n",
        "  x_test = np.asarray(np.load(dirpath + 'x_test.npy', allow_pickle=True)).astype(np.float32)\r\n",
        "  x_valid = np.asarray(np.load(dirpath + 'x_valid.npy', allow_pickle=True)).astype(np.float32)\r\n",
        "  return x_train, x_test, x_valid"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tn-ErKYzhZZx"
      },
      "source": [
        "x_train_unscaled, x_test_unscaled, x_valid_unscaled  = load_datasets()\r\n",
        "print(\"Data loaded. Shapes:\")\r\n",
        "print(x_train_unscaled.shape, x_test_unscaled.shape, x_valid_unscaled.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5EporEztRJ5L"
      },
      "source": [
        "##### Plot distributions of unscaled data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "je6ycNzG3EMw"
      },
      "source": [
        "cols = ['Statistic','F3', 'F4', 'FC3', 'FC4', 'C3', 'Cz', 'C4', 'CP3', 'CP4']\r\n",
        "def data_summary(dataset, label):\r\n",
        "  \"\"\" \r\n",
        "  input:\r\n",
        "    dataset     the three dimensional input (n_samples, n_timepoints, n_features) \r\n",
        "\r\n",
        "    Prints histograms for the 9 features individually\r\n",
        "  returns: \r\n",
        "    summ_df     pd.DataFrame containing summary statistics for the 9 features.\r\n",
        "  \"\"\"\r\n",
        "  data = dataset.reshape((dataset.shape[0] * dataset.shape[1], dataset.shape[2])) # Reshape to 2D (n_samples*n_timepoints, n_features)\r\n",
        "  \r\n",
        "  # Calculate the summary statistics.\r\n",
        "  min   = data.min(axis=0).reshape(1, data.shape[1])                  # Calculate the minimum over the rows for each column.\r\n",
        "  max   = data.max(axis=0).reshape(1, data.shape[1])                  # Then reshape the result to one row and n_cols=n_features, to make it easier to combine later.\r\n",
        "  mean  = data.mean(axis=0).reshape(1, data.shape[1])\r\n",
        "  var   = data.var(axis=0).reshape(1, data.shape[1])\r\n",
        "  q01   = np.quantile(data, 0.01, axis=0).reshape(1, data.shape[1])\r\n",
        "  q99   = np.quantile(data, 0.99, axis=0).reshape(1, data.shape[1])\r\n",
        "\r\n",
        "  names=np.array([['min','max','mean','var','1st percentile', '99th percentile']]).reshape(6,1) # Create a column of names for the summary stats.\r\n",
        "  stats = np.concatenate((min,max,mean,var,q01,q99), axis=0)          # Combine the summary stats in one array\r\n",
        "\r\n",
        "  summ = np.concatenate((names, np.round(stats, 4)), axis=1)          # Combine the summary stats with their names.\r\n",
        "  summ_df = pd.DataFrame(summ, columns=cols)                          # Create a dataframe and supply the channel names as columns.\r\n",
        "\r\n",
        "  # Plot histograms per channel.\r\n",
        "  fig, axes = plt.subplots(3,3, figsize = (9,9))\r\n",
        "  axes=axes.ravel()\r\n",
        "  for i in range(9): # Loop through the channels.\r\n",
        "    axes[i].hist(data[:,i], range= (q01[0,i], q99[0,i]),   density=True)    # Add histogram subplot for the values of that channel.\r\n",
        "    axes[i].title.set_text(cols[i+1])                                       # Add a title with the channel name.\r\n",
        "  fig.suptitle(\"Distribution for each channel (between 1st & 99th percentile) of \" + label + \" dataset\" , size=16)\r\n",
        "  fig.tight_layout(rect=[0, 0.03, 1, 0.95])                                 # Cut the plot space to make space for the global title.\r\n",
        "\r\n",
        "  return summ_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JPhAKBtY04dN"
      },
      "source": [
        "# Plot distributions of each channel.\r\n",
        "unscaled_training_summary = data_summary(x_train_unscaled, \"Training\")\r\n",
        "unscaled_test_summary = data_summary(x_test_unscaled, \"Test\")\r\n",
        "unscaled_valid_summary = data_summary(x_valid_unscaled, \"Validation\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "obsst4fTzSK4"
      },
      "source": [
        "### Band-pass filter"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F6M8UO7FRueQ"
      },
      "source": [
        "##### Create the filters and apply across the whole data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cDMKQR0SzQ47"
      },
      "source": [
        "from scipy import signal\r\n",
        "from copy import deepcopy\r\n",
        "\r\n",
        "low_cut = 0.1\r\n",
        "high_cut = 50.0\r\n",
        "\r\n",
        "bp = signal.butter(10, (low_cut,high_cut), 'bp', fs=500, output='sos') # Create the filter. fs is the sampling rate.\r\n",
        "\r\n",
        "# Create copies of the data\r\n",
        "x_train_filtered = deepcopy(x_train_unscaled)         # After running once in the session, I comment these out because otherwise if you re-run the cell it eats RAM.\r\n",
        "x_test_filtered = deepcopy(x_test_unscaled)\r\n",
        "x_valid_filtered = deepcopy(x_valid_unscaled)\r\n",
        "\r\n",
        "print(x_train_filtered.shape)\r\n",
        "\r\n",
        "x_train_filtered = signal.sosfilt(bp, x_train_filtered, axis=1)\r\n",
        "x_test_filtered = signal.sosfilt(bp, x_test_filtered, axis=1)\r\n",
        "x_valid_filtered = signal.sosfilt(bp, x_valid_filtered, axis=1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0Jttdl9hv6zj"
      },
      "source": [
        "# How many values lay beyond a threshold?\r\n",
        "threshold = 85\r\n",
        "\r\n",
        "outliers_pre_BP = abs(x_train_unscaled) > threshold\r\n",
        "outliers_post_BP = abs(x_train_filtered) > threshold\r\n",
        "                                \r\n",
        "print(\"Before BP filtering, {0:.2f}% of all values lay outside ±{1}\".format((100*np.sum(outliers_pre_BP))/(6665*2500*9), threshold))\r\n",
        "print(\"After BP filtering, {0:.2f}% of all values lay outside ±{1}\".format((100*np.sum(outliers_post_BP))/(6665*2500*9), threshold))\r\n",
        "\r\n",
        " # Number of samples with at least value outside threshold.\r\n",
        "print(\"After BP filtering, there are {0} samples with at least one value outside ±{1}\".format(np.sum(np.max(outliers_post_BP, axis=(1,2))), threshold))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9EQIwyM7tCAp"
      },
      "source": [
        "#### Apply scaling with quantile transformer"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AShIoKEB5_lh"
      },
      "source": [
        "transformer = QuantileTransformer(output_distribution='normal')\r\n",
        "\r\n",
        "# Create copies so the original data remains unaltered.\r\n",
        "x_train = deepcopy(x_train_filtered)\r\n",
        "x_valid = deepcopy(x_valid_filtered)\r\n",
        "x_test = deepcopy(x_test_filtered)\r\n",
        "\r\n",
        "# Get the dimensionality for re-shaping.\r\n",
        "n_samples, n_timepoints, n_features = x_train.shape\r\n",
        "n_samples_valid = x_valid.shape[0]\r\n",
        "n_samples_test = x_test.shape[0]\r\n",
        "\r\n",
        "# Re-shape to 2D for the scaler.\r\n",
        "x_train = x_train.reshape((n_samples*n_timepoints, n_features))        \r\n",
        "x_valid = x_valid.reshape((n_samples_valid*n_timepoints, n_features))\r\n",
        "x_test = x_test.reshape((n_samples_test*n_timepoints, n_features))\r\n",
        "\r\n",
        "# Fit and apply the scaler/transformer to the datasets.\r\n",
        "x_train = transformer.fit_transform(x_train)         \r\n",
        "x_valid = transformer.transform(x_valid)            \r\n",
        "x_test = transformer.transform(x_test)    \r\n",
        "\r\n",
        "# Re-shape to 3D for input to the convolutional autoencoder.\r\n",
        "x_train = x_train.reshape((n_samples, n_timepoints, n_features)) \r\n",
        "x_valid = x_valid.reshape((n_samples_valid, n_timepoints, n_features)) \r\n",
        "x_test = x_test.reshape((n_samples_test, n_timepoints, n_features))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oZQI-k4G6WYu"
      },
      "source": [
        " data_summary(x_train, \"Training Filtered\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dQ7XEKeEo8a5"
      },
      "source": [
        "##### Visualise the Raw & Scaled Signals"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4DEghaTsYrFO"
      },
      "source": [
        "labels=[(('F3',0), ('F4',1)), (('FC3',2), ('FC4',3)), (('C3',4), ('Cz',5), ('C4',6)), (('CP3',7), ('CP4',8))]\r\n",
        "colours= ['darkslateblue', 'orange','lightskyblue','brown','darkgreen','darkgrey','bisque','violet','palegreen']\r\n",
        "\r\n",
        "def plot_signals(sample, title=None):\r\n",
        "  fig, axes = plt.subplots(2,2, figsize = (6,6))\r\n",
        "  axes=axes.ravel()\r\n",
        "  plt.suptitle(\"Signals\" if title is None else title, size=16)\r\n",
        "  count=0\r\n",
        "  for label_group in labels:\r\n",
        "    for label, ind in label_group:\r\n",
        "      axes[count].plot(sample[:,ind], label=label,color=colours[ind], alpha=0.8)\r\n",
        "      axes[count].legend()\r\n",
        "    count+=1\r\n",
        "\r\n",
        "# plot_signals(x_train[101], title=\"Scaled Signals - x_train[0]\")\r\n",
        "plot_signals(x_train_unscaled[101], title=\"Unscaled Signals - x_train[0]\")\r\n",
        "plot_signals(x_train_filtered[101], title=\"Filtered Signals - x_train[0]\")\r\n",
        "plot_signals(x_train[101], title=\"Scaled Signals - x_train[0]\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gG4uVFIx9ZOg"
      },
      "source": [
        "## Utility\r\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B1nJcvib9k7-"
      },
      "source": [
        "#### Experiment Tracker\r\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4-qEtoed9hxZ"
      },
      "source": [
        "report_path = \"/content/drive/MyDrive/ml2-eeg-biometrics/saved_encoders/\"\r\n",
        "class ExperimentTracker():\r\n",
        "\r\n",
        "  def __init__(self):\r\n",
        "    self.last_index = 1\r\n",
        "    try:\r\n",
        "      with open(report_path + \"experiment_tracker.csv\", \"r\") as rfh:\r\n",
        "        reader = csv.reader(rfh, delimiter=\",\", quotechar='\"')\r\n",
        "        # next(reader, None)  # skip the headers\r\n",
        "        data_read = [row for row in reader]\r\n",
        "        print(\"Experiment Tracker found:\", data_read)\r\n",
        "        self.last_index = len(data_read)-2\r\n",
        "        print(\"Last Index: \", self.last_index)\r\n",
        "    except:\r\n",
        "      with open(report_path + \"experiment_tracker.csv\", \"a\") as wfh:\r\n",
        "        print(\"File doesn't exist in path [ \" + report_path +  \" ] === Creating new Experiment Tracker\")\r\n",
        "        self.writer = csv.writer(wfh, delimiter=\",\")\r\n",
        "        # writer.writerow([\"your\", \"header\", \"foo\"])  # write header\r\n",
        "        header = [\"Experiment N°\", \"Model Name\", \"Kernel Sizes\", \"Filters\", \"Pool Size\", \"Size of Encoding\", \"Loss\"]\r\n",
        "        self.writer.writerow(header)\r\n",
        "    with open(report_path + \"experiment_tracker.csv\", \"a\") as wfh:\r\n",
        "      self.writer = csv.writer(wfh, delimiter=\",\")\r\n",
        "\r\n",
        "  def report_result(self, model_name, kernel_sizes, filters, pool_size, encoding_size, train_loss, val_loss, train_acc, val_acc):\r\n",
        "    with open(report_path + \"experiment_tracker.csv\", \"a\") as wfh:\r\n",
        "      self.last_index += 1\r\n",
        "      writer = csv.writer(wfh, delimiter=\",\")\r\n",
        "      row = [str(self.last_index), str(model_name), str(kernel_sizes), str(filters), str(pool_size), str(encoding_size), str(train_loss[0]), str(val_loss[0]), str(train_acc[0]), str(val_acc[0])]\r\n",
        "      writer.writerow(row)\r\n",
        "      print(\"Row saved:\", row)\r\n",
        "\r\n",
        "  def get_last_index(self):\r\n",
        "    return self.last_index\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I7PWDBiuWVX8"
      },
      "source": [
        "## Autoencoders"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UIUsOKCEyxv1"
      },
      "source": [
        "#### Convolutional Autoencoder"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ltaV2CiYXS5J"
      },
      "source": [
        "##### Define a convolutional autoencoder."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N1xX-IPc1GXH"
      },
      "source": [
        "This object-orientd implementation of the Convolutional Autoencoder allows to build an extra candidate architecture to be compared with the current architecture. In additions it includes several methods (that wrap native methdos from the class Model) to plot, evaluate and save a model easily while performing an experiment."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WjpduhZRc2TF"
      },
      "source": [
        "path_models = \"/content/drive/MyDrive/ml2-eeg-biometrics/saved_models/\"\r\n",
        "path_encoders = \"/content/drive/MyDrive/ml2-eeg-biometrics/saved_encoders/\"\r\n",
        "\r\n",
        "class ConvAutoencoder():\r\n",
        "\r\n",
        "  def __init__(self, loss, optimizer, shape, parameters, tracker, candidate=False):\r\n",
        "      self.input_shape = shape   \r\n",
        "      self.encoded_layer = None\r\n",
        "      self.history = None\r\n",
        "      self.tracker = tracker\r\n",
        "      self.parameters = parameters\r\n",
        "      print(\"Candidate build:\", candidate)    \r\n",
        "      self.autoencoder_model = self.build_candidate_model() if candidate else self.build_model(parameters['kernel_size'], parameters['filters'], parameters['pool_size'])\r\n",
        "      self.encoder = Model(inputs=self.autoencoder_model.input, outputs=self.encoded_layer) # save the encoder for later extraction\r\n",
        "      self.autoencoder_model.compile(loss=loss, optimizer=optimizer)\r\n",
        "      self.autoencoder_model.summary()\r\n",
        "  \r\n",
        "\r\n",
        "  ''' Builds the architecture of the model'''   \r\n",
        "  def build_model(self, kernel_sizes, filters, pool_sizes):\r\n",
        "      input_layer = Input(shape=self.input_shape)\r\n",
        "\r\n",
        "      # encoder\r\n",
        "      e1 = Conv1D(filters=filters[0], kernel_size=kernel_sizes[0], activation='relu', padding='same')(input_layer)\r\n",
        "      e2 = MaxPooling1D(pool_size=pool_sizes[0])(e1)\r\n",
        "      e3 = Conv1D(filters=filters[1], kernel_size=kernel_sizes[1], activation = 'relu', padding='same')(e2)\r\n",
        "      e4 = MaxPooling1D(pool_size=pool_sizes[1])(e3)\r\n",
        "      e5 = Conv1D(filters=filters[2], kernel_size=kernel_sizes[2], activation = 'relu', padding='same')(e4)\r\n",
        "      e6 = MaxPooling1D(pool_size=pool_sizes[2])(e5)\r\n",
        "      encoded = Flatten()(e6)\r\n",
        "      self.encoded_layer = encoded\r\n",
        "      encoding_size = encoded.type_spec.shape[1]\r\n",
        "      print(\"Encoder size: \", encoding_size)\r\n",
        "\r\n",
        "      # decoder\r\n",
        "      reshape_size = int(encoding_size/2)\r\n",
        "      d0 = Reshape((reshape_size, 2))(encoded)\r\n",
        "      d1 = UpSampling1D(size=pool_sizes[2])(d0)\r\n",
        "      d2 = Conv1DTranspose(filters=filters[1], kernel_size=kernel_sizes[2], activation='relu', padding='same')(d1)    \r\n",
        "      d3 = UpSampling1D(size=pool_sizes[1])(d2)\r\n",
        "      d4 = Conv1DTranspose(filters=filters[2], kernel_size=kernel_sizes[1], activation='relu', padding='same')(d3)\r\n",
        "      d5 = UpSampling1D(size=pool_sizes[0])(d4)\r\n",
        "      decoded = Conv1DTranspose(filters=9, kernel_size=kernel_sizes[0], activation='linear', padding='same')(d5)\r\n",
        "\r\n",
        "      model = Model(inputs=input_layer, outputs=decoded)\r\n",
        "      model.output_shape\r\n",
        "      return model\r\n",
        "\r\n",
        "\r\n",
        "  ''' Builds the architecture of a candidate model, to be used in comparison with the current model'''   \r\n",
        "  def build_candidate_model(self):\r\n",
        "      print(\"No candidate build present, implement one first\")\r\n",
        "      return None\r\n",
        "\r\n",
        "\r\n",
        "  ''' Trains the model on the full dataset, requires validation split'''                     \r\n",
        "  def train_model(self, x_train, x_val, epochs, batch_size=20):\r\n",
        "      early_stopping = EarlyStopping(monitor='loss',\r\n",
        "                                      min_delta=0,\r\n",
        "                                      patience=5,\r\n",
        "                                      verbose=1, \r\n",
        "                                      mode='auto')\r\n",
        "      self.history = self.autoencoder_model.fit(x_train, x_train,\r\n",
        "                                            batch_size=batch_size,\r\n",
        "                                            epochs=epochs,\r\n",
        "                                            validation_data=(x_val, x_val),\r\n",
        "                                            callbacks=[early_stopping])\r\n",
        "      return self.history\r\n",
        "\r\n",
        "\r\n",
        "  ''' Trains the model on a single channel, no validation split'''   \r\n",
        "  def train_subset(self, x_train, epochs, batch_size=20):\r\n",
        "      history = self.autoencoder_model.fit(x_train[:,:,0], x_train[:,:,0], epochs=epochs,\r\n",
        "                  batch_size=batch_size,\r\n",
        "                  # validation_split=0.05\r\n",
        "                  )\r\n",
        "      \r\n",
        "\r\n",
        "  ''' Plots the history.\r\n",
        "  input:\r\n",
        "    history: if present, plot the history given as parameter\r\n",
        "  output:\r\n",
        "    None\r\n",
        "  '''       \r\n",
        "  def plot_history(self, history=None):\r\n",
        "      history = history if history is not None else self.history\r\n",
        "      plt.plot(history.history['loss'])\r\n",
        "      plt.plot(history.history['val_loss'])\r\n",
        "      plt.title('Model loss')\r\n",
        "      plt.ylabel('Loss')\r\n",
        "      plt.xlabel('Epoch')\r\n",
        "      plt.legend(['Train', 'Test'], loc='upper left')\r\n",
        "      plt.show()\r\n",
        "      \r\n",
        "\r\n",
        "  ''' Return the predictions for the trained data.\r\n",
        "   input:\r\n",
        "    X: the training data\r\n",
        "   output:\r\n",
        "    preds: predictions for the given training data\r\n",
        "  '''       \r\n",
        "  def eval_model(self, X):\r\n",
        "      preds = self.autoencoder_model.predict(X)\r\n",
        "      return preds\r\n",
        "\r\n",
        "\r\n",
        "  ''' Save all model information, including weights, in h5 format\r\n",
        "    input:\r\n",
        "      path: full path on drive, for example /content/drive/MyDrive/ml2-eeg-biometrics/saved_models/\r\n",
        "    output:\r\n",
        "      None\r\n",
        "  '''    \r\n",
        "  def save_model(self, path, label=''):\r\n",
        "      timestamp = pd.Timestamp.now()\r\n",
        "      model_name = label + \"_\" + str(timestamp) + \".h5\"\r\n",
        "      saved_model = self.autoencoder_model.save(path + model_name)\r\n",
        "      print(\"Model saved at \" + str(timestamp) + \" in path \" + path)\r\n",
        "\r\n",
        "  \r\n",
        "  ''' Writes a line with metadata in the Experiment Tracker, then save the encoder part of the model, \r\n",
        "      including weights, in h5 format. Finally, saves the encodings in .npy format.\r\n",
        "    input:\r\n",
        "      path: full path on drive, for example /content/drive/MyDrive/ml2-eeg-biometrics/saved_models/\r\n",
        "      label: name of the file\r\n",
        "      train_pred: predictions for the training set on the encoder\r\n",
        "      valid_pred: predictions on the validation set on the encoder\r\n",
        "      loss: loss of the model\r\n",
        "    output:\r\n",
        "      None\r\n",
        "  '''    \r\n",
        "  def report_experiment_result(self, path, label, train_pred, valid_pred):\r\n",
        "      timestamp = pd.Timestamp.now()\r\n",
        "      model_name = label + \"_\" + str(timestamp) \r\n",
        "      train_loss, val_loss = self.retrieve_losses()    \r\n",
        "      train_acc, val_acc = self.retrieve_accuracies()\r\n",
        "      self.tracker.report_result(model_name + \".h5\", self.parameters[\"kernel_size\"], self.parameters[\"filters\"], self.parameters[\"pool_size\"],  self.encoder.output_shape[1], train_loss, val_loss,  train_acc, val_acc)\r\n",
        "      \r\n",
        "      index = self.tracker.get_last_index()\r\n",
        "      saved_encoder = self.encoder.save(path + str(index) + \"_\" +  model_name + \".h5\")\r\n",
        "      print(\"Encoder saved at \" + str(timestamp) + \" in path \" + path)\r\n",
        "      np.save(path + \"encodings/\" + str(index) + \"_train_encoding_\" + model_name + \".npy\", train_pred)\r\n",
        "      np.save(path + \"encodings/\" + str(index) + \"_valid_encoding_\" + model_name + \".npy\", valid_pred)\r\n",
        "\r\n",
        "  ''' Retrieve the history of the model.\r\n",
        "  input:\r\n",
        "    None\r\n",
        "  output:\r\n",
        "    history: history saved in the object model\r\n",
        "  '''\r\n",
        "  def retrieve_history(self):\r\n",
        "      return self.history\r\n",
        "\r\n",
        "  ''' Retrieve the losses of the model.\r\n",
        "  input:\r\n",
        "    None\r\n",
        "  output:\r\n",
        "    train_loss: the loss on the training data\r\n",
        "    val_loss: the loss on the validation data\r\n",
        "  '''\r\n",
        "  def retrieve_losses(self):\r\n",
        "    if 'loss' in self.history.history and 'val_loss' in self.history.history :\r\n",
        "      train_loss = self.history.history['loss']\r\n",
        "      val_loss   = self.history.history['val_loss']\r\n",
        "      return train_loss, val_loss\r\n",
        "    else:\r\n",
        "      return [''], ['']\r\n",
        "\r\n",
        "  ''' Retrieve the accuracies of the model.\r\n",
        "  input:\r\n",
        "    None\r\n",
        "  output:\r\n",
        "    train_acc: the accuracy on the training data\r\n",
        "    val_acc: the accuracy on the validation data\r\n",
        "  '''\r\n",
        "  def retrieve_accuracies(self):\r\n",
        "    if 'acc' in self.history.history and 'val_acc' in self.history.history :\r\n",
        "      train_acc  = self.history.history['acc']\r\n",
        "      val_acc    = self.history.history['val_acc']\r\n",
        "      return train_acc, val_acc\r\n",
        "    else:\r\n",
        "      return [''], ['']\r\n",
        "\r\n",
        "  ''' Retrieve the encoder part of the model.\r\n",
        "  input:\r\n",
        "    index: the index of the encoded layer\r\n",
        "  output:\r\n",
        "    encoder: a model built with the input layer as input and the encoded layer as output\r\n",
        "  '''\r\n",
        "  def retrieve_encoder(self, index):\r\n",
        "      return self.encoder\r\n",
        "  \r\n",
        "\r\n",
        "  ''' Assess the model on a subset of data\r\n",
        "  input:\r\n",
        "    X: the training data\r\n",
        "  output:\r\n",
        "    None\r\n",
        "  '''\r\n",
        "  def assess_subset(self, X):\r\n",
        "    x_pred = self.eval_model(X)\r\n",
        "    for i in range(0,len(x_1),200):\r\n",
        "      plt.plot(x_1[i,:,0], label='actual')\r\n",
        "      plt.plot(x_pred[i,:,0], label = 'predicted')\r\n",
        "      plt.show()\r\n",
        "\r\n",
        "\r\n",
        "  ''' Plot prediction for the model.\r\n",
        "  input:\r\n",
        "    X: the training data\r\n",
        "  output:\r\n",
        "    None\r\n",
        "  '''\r\n",
        "  def plot_prediction(self, X):\r\n",
        "    x_pred = candidate_model.eval_model(X)\r\n",
        "    plt.figure(figsize=(14,8))\r\n",
        "    plt.plot(X[1039,:,0], label='actual', alpha=0.7)\r\n",
        "    plt.plot(x_pred[1039,:,0], label= 'predicted',alpha=0.7)\r\n",
        "    # plt.ylim(0.475, 0.495)\r\n",
        "    plt.legend()\r\n",
        "\r\n",
        "  ''' Plot the distribution of the loss across the dataset\r\n",
        "  input:\r\n",
        "    X_pred: prediction on the training data (obtainable from eval_model())\r\n",
        "    X_train: the training data\r\n",
        "  output:\r\n",
        "    None\r\n",
        "  '''\r\n",
        "  def plot_loss_distribution(self, X_pred, x_train):\r\n",
        "    x_train_reshaped = x_train.reshape(x_train.shape[0]*x_train.shape[1], x_train.shape[2])\r\n",
        "    loss_mae = np.mean(np.abs(X_pred-x_train_reshaped), axis = 1)\r\n",
        "    plt.figure(figsize=(16,9), dpi=80)\r\n",
        "    plt.title('Loss Distribution', fontsize=16)\r\n",
        "    sns.distplot(loss_mae, bins = 20, kde= True, color = 'blue');\r\n",
        "    plt.xlim([0.0,.5])\r\n",
        "    \r\n",
        "\r\n",
        "  \"\"\" Function to plot predictions vs. the actuals for one sample.\r\n",
        "  input:\r\n",
        "    actuals   3D array (n_samplesx2500x9) - Original scaled signals.\r\n",
        "    pred      the predicted values corresponding to the actuals.\r\n",
        "    ind     The row number of the sample (2500x9) that you want to compare.\r\n",
        "    rescale   If set to true then the data is first converted back to the original scale for comparison.\r\n",
        "  \r\n",
        "  returns:\r\n",
        "    nothing\r\n",
        "\r\n",
        "  prints plots.\r\n",
        "  \"\"\"\r\n",
        "  def evaluate_prediction(self, actuals, pred, ind, rescale=False):\r\n",
        "    cols = ['F3', 'F4', 'FC3', 'FC4', 'C3', 'Cz', 'C4', 'CP3', 'CP4']\r\n",
        "\r\n",
        "    if rescale:\r\n",
        "      sample_actual = scaler.inverse_transform(actuals[ind]) # Rescale to the original scale.\r\n",
        "      sample_pred = scaler.inverse_transform(pred[ind]) \r\n",
        "    else:\r\n",
        "      sample_actual = actuals[ind]    # Get the relevant sample.\r\n",
        "      sample_pred = pred[ind]\r\n",
        "\r\n",
        "    mae_by_channel = np.mean(np.abs(sample_pred - sample_actual), axis=0) # Get the Mean Absolute Error for each channel for this sample\r\n",
        "    sample_mae = np.mean(mae_by_channel) # Get the total MAE for the sample by taking the average across the 9 channels\r\n",
        "    print(\"Sample\", ind, \"\\n   Total Mean Absolute Error:\", round(sample_mae, 8))\r\n",
        "    print(\"Mean Absolute Error by Channel:\")\r\n",
        "    for col, error in zip(cols, mae_by_channel):\r\n",
        "      print(col, \": \", round(error,8)) \r\n",
        "\r\n",
        "    fig, axes = plt.subplots(3,3, figsize=(9,9))\r\n",
        "    axes=axes.ravel()\r\n",
        "\r\n",
        "    for i in range(9):\r\n",
        "      axes[i].plot(sample_actual[:,i], label= \"Actual\")\r\n",
        "      axes[i].plot(sample_pred[:,i], label=\"Predicted\")\r\n",
        "      axes[i].title.set_text(cols[i] + str(round(mae_by_channel[i], 3)))\r\n",
        "    \r\n",
        "    plt.legend()\r\n",
        "    fig.suptitle(\"Predictions vs. Actuals - Sample \" + str(ind),size=16)\r\n",
        "    fig.tight_layout(rect=[0, 0.03, 1, 0.95])\r\n",
        "\r\n",
        "  ''' Plot loss distribution of the test set for each channel.\r\n",
        "  input:\r\n",
        "    X_test: the test data \r\n",
        "  output:\r\n",
        "    None\r\n",
        "  '''\r\n",
        "  def plot_loss_test(self, x_test):\r\n",
        "    X_pred = self.eval_model(x_test)\r\n",
        "    X_pred = X_pred.reshape(X_pred.shape[0]*X_pred.shape[1], X_pred.shape[2])\r\n",
        "\r\n",
        "    x_test_reshaped = x_test.reshape(x_test.shape[0]*x_test.shape[1], x_test.shape[2])\r\n",
        "    fig, axes = plt.subplots(9,1, figsize=(18,9))\r\n",
        "    # Plot the loss distribution for each channel individually\r\n",
        "    for i in range(x_test_reshaped.shape[1]):\r\n",
        "      loss_mae = np.abs(X_pred[:,i]-x_test_reshaped[:,i])\r\n",
        "      sns.distplot(loss_mae, bins = 100, kde= True, color = 'blue', ax=axes[i]);\r\n",
        "      axes[i].axis(xmin=0.0,xmax=0.2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SUwCCZ0ZXxdR"
      },
      "source": [
        "##### Fit the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YbFLfZEL_Gdh"
      },
      "source": [
        "# create the autoencoder model with the current architecture\r\n",
        "tracker = ExperimentTracker()\r\n",
        "optimizer = Adam(lr=0.001)\r\n",
        "shape = (x_train.shape[1], x_train.shape[2])\r\n",
        "parameters = {\r\n",
        "          \"kernel_size\": [9, 5, 3],\r\n",
        "          \"filters\": [18, 6, 2],\r\n",
        "          \"pool_size\": [2, 5, 2]\r\n",
        "      }\r\n",
        "      \r\n",
        "model = ConvAutoencoder(loss='mae', optimizer=optimizer, shape=shape, parameters=parameters, tracker=tracker)\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fzzXJTvNyAn0"
      },
      "source": [
        "# create the autoencoder model with the candidate architecture\r\n",
        "#optimizer = Adam(lr=0.001)\r\n",
        "#shape = (x_train.shape[1], x_train.shape[2])\r\n",
        "#candidate_model = ConvAutoencoder(loss='mae', optimizer=optimizer, shape=shape, candidate=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZDgOztHOHazh"
      },
      "source": [
        "# Train both current and candidate models\r\n",
        "\r\n",
        "# fit the model to the data\r\n",
        "nb_epochs = 1\r\n",
        "batch_size = 20\r\n",
        "x_1 = x_train[:,:,:]\r\n",
        "## Train the models on the full dataset\r\n",
        "print('Training Current model for ' + str(nb_epochs) + \" epochs and batch size of \" + str(batch_size))\r\n",
        "model.train_model(x_1, x_valid, nb_epochs, batch_size)\r\n",
        "print('============================================================')\r\n",
        "#print('Training Candidate model for ' + str(nb_epochs) + \" epochs and batch size of \" + str(batch_size))\r\n",
        "#candidate_model.train_model(x_1, x_valid, nb_epochs, batch_size)\r\n",
        "#print('============================================================')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "64K4Tw5CcXZs"
      },
      "source": [
        "\r\n",
        "# encoder.summary()\r\n",
        "\r\n",
        "train_pred = model.encoder.predict(x_train)\r\n",
        "valid_pred = model.encoder.predict(x_valid)\r\n",
        "\r\n",
        "model.report_experiment_result(path_encoders, \"encoder_\" + str(nb_epochs) +\"epochs\" + \"_batchsize\" + str(batch_size), train_pred, valid_pred)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LKottLNgXpw6"
      },
      "source": [
        "model.plot_history()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YSGDjFwgX8pZ"
      },
      "source": [
        "##### Plot the reconstruction for individual samples.\r\n",
        "Temporary code for assessing the models on the subset of data (one channel and only 20% of timepoints)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8Hm2DWhiSmeV"
      },
      "source": [
        "model.assess_subset(x_1)\r\n",
        "#candidate_model.assess_subset(x_1)\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yefRYH1sc5LY"
      },
      "source": [
        "Plot prediction for the whole dataset"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "abbHkzSliPCa"
      },
      "source": [
        "model.plot_prediction(x_1)\r\n",
        "#candidate_model.plot_prediction(x_1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kCfIG-Na0Rim"
      },
      "source": [
        "\r\n",
        "## Evaluate Results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x0AvXQtdVLPT"
      },
      "source": [
        "# Get the predicted values for the training set.\r\n",
        "X_pred_3D = model.eval_model(x_train)\r\n",
        "# X_pred_3D = candidate_model.eval_model(x_train_f)\r\n",
        "X_pred = X_pred_3D.reshape(X_pred_3D.shape[0]*X_pred_3D.shape[1], X_pred_3D.shape[2])\r\n",
        "# X_pred = pd.DataFrame(X_pred, columns=train.columns)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IIl53WipVnBA"
      },
      "source": [
        "##### Plot distribution of the loss"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OERj4N100Maz"
      },
      "source": [
        "# Plot the distribution of the loss\r\n",
        "model.plot_loss_distribution(X_pred, x_train)\r\n",
        "# candidate_model.plot_loss_distribution(X_pred, x_train)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5uQFT7_1Vtt1"
      },
      "source": [
        "##### Evaluate total re-construction for one sample"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UKYRpNQjWF-u"
      },
      "source": [
        "model.evaluate_prediction(x_train, X_pred_3D, ind=1, rescale=False)\r\n",
        "# candidate_model.evaluate_prediction(x_train, X_pred_3D, ind=1, rescale=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-ehGjGxUV4Vf"
      },
      "source": [
        "##### Plot the loss distribution of the test set (code from tutorial, needs to be edited)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sk2FB1JADjIE"
      },
      "source": [
        "# plot the loss distribution of the test set\r\n",
        "model.plot_loss_test(x_test)\r\n",
        "# candidate_model.plot_loss_test(x_test)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XVIfqljx0lwt"
      },
      "source": [
        "model.save_model(path_models, \"model_20epochs\") # Get current timestamp and model number.\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Wbha7pImw6J2"
      },
      "source": [
        "#Grid Search Optimizer"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XDYqnExxw98w"
      },
      "source": [
        "#### This algorithm tries 9 possibile combination of hyperparameters and saves the encoder model and the encoding vectors as files in the drive folder"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R3pj58nwxO3T"
      },
      "source": [
        "\r\n",
        "\r\n",
        "class GridSeachExperimentalRunner():\r\n",
        "\r\n",
        "  def __init__(self, x_train, x_valid):\r\n",
        "    self.x_train = x_train\r\n",
        "    self.x_valid = x_valid\r\n",
        "    self.kernel_sizes = [[9, 5, 3], [11, 7, 5], [5, 5, 5]]\r\n",
        "    self.filters = [[18, 6, 2], [27, 18, 2], [27, 9, 2]]\r\n",
        "    self.pool_sizes = [[2, 5, 2], [2, 5, 1], [2, 2, 1]]\r\n",
        "    self.nb_epochs = [10, 25, 50]\r\n",
        "    self.batch_sizes = [20, 50, 100]\r\n",
        "    self.models = []\r\n",
        "\r\n",
        "  def build_models(self):\r\n",
        "    # create the autoencoder model with the current architecture\r\n",
        "    tracker = ExperimentTracker()\r\n",
        "    optimizer = Adam(lr=0.001)\r\n",
        "    shape = (x_train.shape[1], x_train.shape[2])\r\n",
        "    for kernel_size in self.kernel_sizes:\r\n",
        "      for filter in self.filters:\r\n",
        "        for pool_size in self.pool_sizes:\r\n",
        "          parameters = {\r\n",
        "                    \"kernel_size\": kernel_size,\r\n",
        "                    \"filters\": filter,\r\n",
        "                    \"pool_size\": pool_size\r\n",
        "                }   \r\n",
        "          model = ConvAutoencoder(loss='mae', optimizer=optimizer, shape=shape, parameters=parameters, tracker=tracker) \r\n",
        "          self.models.append(model)\r\n",
        "    print(\"Candidate models: \", len(self.models))\r\n",
        "\r\n",
        "  def summarize_model(self, model, x_train, x_valid, nb_epochs, batch_size):\r\n",
        "    train_pred = model.encoder.predict(x_train)\r\n",
        "    valid_pred = model.encoder.predict(x_valid)\r\n",
        "    model.report_experiment_result(path_encoders, \"encoder_\" + str(nb_epochs) +\"epochs\" + \"_batchsize\" + str(batch_size), train_pred, valid_pred) \r\n",
        "\r\n",
        "  def run_experiment(self):\r\n",
        "    print('Training candidate models for ' + str(self.nb_epochs[2]) + \" epochs and batch size of \" + str(self.batch_sizes[0]))\r\n",
        "\r\n",
        "    for model in self.models:\r\n",
        "      model.train_model(x_train, x_valid, self.nb_epochs[2], self.batch_sizes[0])\r\n",
        "      self.summarize_model(model, x_train, x_valid, self.nb_epochs[2], self.batch_sizes[0])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6Sq1hZvz4GVc"
      },
      "source": [
        "grid_search = GridSeachOptimizer(x_train, x_valid)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lJHNnPMX0qz0"
      },
      "source": [
        "# WARNING: this code will produce 27 models and relative summary output\r\n",
        "\r\n",
        "grid_search.build_models()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DT1cbRvN24Yl"
      },
      "source": [
        "# WARNING: this code will train 27 models for 50 epochs and batch size of 20 (default values)\r\n",
        "grid_search.run_experiment()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EpfXOsm8mC_k"
      },
      "source": [
        "# loading the saved encoder\r\n",
        "loaded_encoder = tf.keras.models.load_model(path_encoders + '17_encoder_25epochs_batchsize20_2021-01-21 22:41:29.870146.h5')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gyP1ZrDAnOr2"
      },
      "source": [
        "loaded_encoder.summary()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4wH8_MFfnTg1"
      },
      "source": [
        "# predict on test data and then save\r\n",
        "test_pred = loaded_encoder.predict(x_test)\r\n",
        "np.save(path_encoders + \"test_encodings/\" + \"17_test_encoding.npy\", test_pred)\r\n",
        "\r\n"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}