{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification of cardiomyocites based on their action potentials\n",
    "In this lab experience we will be looking at the problem of classifying cardiac cells by looking at\n",
    "their _Action Potentials_ (APs). For the purpose of this task we will be using synthetically\n",
    "generated APs following the models in [(Nygren et. al, 1998)](https://www.ahajournals.org/doi/abs/10.1161/01.res.82.1.63) and [(O'Hara et. al, 2011)](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002061).\n",
    "\n",
    "### Objectives\n",
    "By the end of this session you should be able to use nearest-neighbor and simple neural network classifiers for two-class classification problems.\n",
    "\n",
    "* _Note: Everything below this line should be truned in as part of your lab report._\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============================================================================\n",
    "# import modules\n",
    "# ============================================================================\n",
    "# Note that this part of the code needs to be run prior to any other code cell\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io\n",
    "import random as rnd\n",
    "import torch\n",
    "\n",
    "# inline plots\n",
    "%matplotlib inline\n",
    "\n",
    "# set random number generator seed\n",
    "rnd.seed(10)\n",
    "\n",
    "# mount GDrive\n",
    "from google.colab import drive\n",
    "drive.mount('gdrive/')\n",
    "\n",
    "# ============================================================================\n",
    "# data loading\n",
    "# ============================================================================\n",
    "\n",
    "# PLEASE, SPECIFY YOUR PATH TO THE DATA HERE\n",
    "datapath = 'gdrive/My Drive/bmdslab/lab-02/'\n",
    "\n",
    "# get the list of files to process\n",
    "matfile = '/'.join((datapath,'Adult_samples.mat'))\n",
    "\n",
    "# get data as dictionary\n",
    "adata = scipy.io.loadmat(matfile)\n",
    "\n",
    "# see the keys of the dictionary 'Vatrial' and 'Vventricular'\n",
    "# adata['Vatrial'] contains action potentials of cardiac cells of atrial type\n",
    "# adata['Vventricular'] contains action potentials of cardiac cells of ventricular type\n",
    "adata.keys()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# look at the data dimensions\n",
    "adata['Vatrial'].shape\n",
    "\n",
    "# let's see how the data looks like\n",
    "plt.plot(adata['Vatrial'][:,0])\n",
    "plt.plot(adata['Vventricular'][:,0])\n",
    "plt.legend(('atrial','ventricular'))\n",
    "plt.ylabel('Voltage [mV]')\n",
    "plt.xlabel('sample')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Task 1. Pre-processing and manual feature extraction\n",
    "In this first part we will be normalizing the data and extracting hand-crafted features that will be later\n",
    "used for the classification task. Load the dataset **Adult_samples.mat** that contains examples of\n",
    "adult atrial and ventricular action potentials generated using the models in [(Nygren et al., 1998)](https://www.ahajournals.org/doi/abs/10.1161/01.res.82.1.63) and [(O'Hara et al., 2011)](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002061).\n",
    "For each of the two classes there are $1000$ samples generated with a sampling rate of $f_s=500$ Hz.\n",
    "\n",
    "* _(15 points) Data preparation and normalization._\n",
    "Split the data between test and training sets by\n",
    "randomly selecting $10\\%$ of the points as your training set. Make sure the two classes are well represented\n",
    "in the training set (e.g., use the same number for both). Normalize the data so that each AP has zero resting potential and unit maximum amplitude. Create an array of corresponding labels for the data points. For ventricular type use the class label $+1$ and $-1$ for\n",
    "atrial type.\n",
    "Make two plots displaying your normalized training data for each of the classes. Use time units for the horizontal axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* _(20 points) Hand-crafted features._\n",
    "The _Action Potential Duration_ (APD) at $x\\%$ is defined as the time\n",
    "it takes to reduce the maximum amplitude of the AP to $x\\%$ of its value. Write a function that computes APD at a given\n",
    "percentage $x\\in[0,1]$. Compute also the _Average of the Action Potential_ (AAP) and build two-dimensional features\n",
    "by concatenating APD@$0.5$ and AAP. Make a scatter plot of the training data using these two features. Use\n",
    "different colors and/or markers to represent each class.\n",
    "\n",
    "    * Based on your scatter plot, is the training data using the above features linearly separable? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============================================================================\n",
    "# compute hand-crafted features such as APA and APDx\n",
    "# ============================================================================\n",
    "\n",
    "def apdx(X,alpha=0.5,fs=1.0):\n",
    "    \"\"\"\n",
    "    This function computes the APDx feature defined as the time it takes to\n",
    "    reduce the action potential to x% of its maximum value.\n",
    "    Use:\n",
    "        APD = apdx(X,alpha,fs)\n",
    "        \n",
    "    Inputs:\n",
    "        X : data points (each column is a data point)\n",
    "        alpha : [0,1] percentage of maximum amplitude\n",
    "        fs : sampling frequency (to map samples to time)\n",
    "        \n",
    "    Output:\n",
    "        APD feature\n",
    "    \"\"\"\n",
    "\n",
    "    # write your function code here\n",
    "    # ...\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2. Classification\n",
    "\n",
    "* _(15 points) Nearest-neighbor classifier._\n",
    "Implement a $1$NN classifier using the Euclidean distance. A $1$NN classifier works as follows: Given your training dataset $\\mathcal{D} = \\big\\{(\\boldsymbol x_i,y_i)\\big\\}_{i=1}^{N}$,where $N$ is the number of training samples, $\\boldsymbol x_i\\in\\mathbb{R}^D$ is a feature vector and $y_i\\in\\{-1,1\\}$ its associated label, and a novel sample $\\boldsymbol x$, the $1$NN classifier assigns to $\\boldsymbol x$ the same label as its closest point in the training set. That is, the estimated label $\\widehat y$ of $\\boldsymbol x$ is such that:\n",
    "\n",
    "\\begin{equation}\n",
    "\t\\widehat y(\\boldsymbol x) = y_{k^*},\\quad k^* = \\underset{i\\in\\{1,\\ldots,N\\}}{\\arg\\min}\\, \\lVert \\boldsymbol x_i - \\boldsymbol x \\rVert_2.\n",
    "\\end{equation}\n",
    "\n",
    "Compute and display the classification accuracy over the test set using the handcrafted training features of **Task 1**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============================================================================\n",
    "# Nearest-neighbor classifier\n",
    "# ============================================================================\n",
    "\n",
    "def nn_classify(Xtest,Xtrain,ytrain):\n",
    "    \"\"\"\n",
    "    This function implements a nearest neighbor classifier.\n",
    "    Use:\n",
    "        yhat = nn_classify(Xtest,Xtrain,ytrain)\n",
    "        \n",
    "    Inputs:\n",
    "        Xtest  : test data points (each column is a data point)\n",
    "        Xtrain : training data points\n",
    "        ytrain : associated labels to the training data points\n",
    "        \n",
    "    Output:\n",
    "        yhat : estimated labels for the test data.\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    # implement your NN classifier below\n",
    "    # ...\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* _(40 points) Two-layer Neural Network._\n",
    " Implement a two-layer neural network classifier of the form:\n",
    "\n",
    "   $$\\widehat y = \\operatorname{sign}\\big(f_{\\theta}(\\phi(\\boldsymbol x))\\big),\\quad f_{\\theta}(\\phi(\\boldsymbol x)) = \\boldsymbol w^T\\phi(\\boldsymbol x) + b,\\quad \\theta = \\left[\\begin{array}{c}\\boldsymbol w\\\\b\\end{array}\\right],$$\n",
    "  \n",
    " where $f_\\theta(\\cdot)$ is a _linear prediction function_ (i.e., classification layer) parametrized by $\\theta = [\\boldsymbol w^T, b]^T$. The feature extraction part of the network consists of a linear layer followed by a ReLu (rectified linear unit) non-linearity:\n",
    " \n",
    " $$ \\phi(x) = \\operatorname{ReLu}\\big(\\boldsymbol W_1 \\boldsymbol x + \\boldsymbol b_1\\big),\\quad\\operatorname{ReLu}(x)=\\begin{cases}x &x>0\\\\0&\\textrm{else}\\end{cases}.$$\n",
    " \n",
    " In order to find the network's parameters $\\Theta=\\{\\boldsymbol W_1, \\boldsymbol b_1, \\boldsymbol w, b\\}$ minimize the following _regularized empirical risk_ using PyTorch:\n",
    "\n",
    "   $$\\min_{\\Theta}\\, \\underbrace{\\frac{1}{N} \\sum_{i=1}^N L\\big(f_\\theta(\\phi(\\boldsymbol x_i)),y_i\\big) + \\lambda\\big(\\frac{1}{w}\\lVert \\boldsymbol w\\rVert^2 + \\frac{1}{W}\\lVert\\boldsymbol W_1\\rVert^2\\big)}_{C(\\Theta)}$$\n",
    "    \n",
    " where the loss function $L(f, y) = \\lVert y-f\\rVert_2^2)$ is the quadratic (square) loss, and $w,W$, are the number of elements in $\\boldsymbol w$ and $\\boldsymbol W$, respectively. \n",
    "\n",
    "    * Define a network model in PyTorch according to the definition above. For that purpose you can use `torch.nn.Sequential`. Follow this [example](https://pytorch.org/tutorials/beginner/pytorch_with_examples.html#pytorch-nn) to learn how to use it.\n",
    "    \n",
    "    * Run a gradient descent algorithm to minimize the cost function using $\\lambda =1$. Carefully choose the stepsize and number of iterations until you see the method converges (i.e., the cost function gets to a \"plateau\").\n",
    "    \n",
    "    * Make a scatter plot of the learned features (i.e., prior to classificaiton layer) by your network model. Has your model learned features that are linearly separable? Display in the scatter plot the decision boundary that you have learned. Compute the classification accuracy over the test set.\n",
    "    \n",
    "    * Plot the weights of the learned linear layer. What has your network learned?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}