unit04_model_sel/lab_neural_partial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab:  Model Order Selection for Neural Data\n",
    "\n",
    "Machine learning is a key tool for neuroscientists to understand how sensory and motor signals are encoded in the brain.  In addition to improving our scientific understanding of neural phenomena, understanding neural encoding is critical for brain machine interfaces.  In this lab, you will use model selection for performing some simple analysis on real neural signals.  \n",
    "\n",
    "Before doing this lab, you should review the ideas in the [polynomial model selection demo](./polyfit.ipynb).  In addition to the concepts in that demo, you will learn to:\n",
    "* Represent neural time-series data in arrays\n",
    "* Load data from a pickle file\n",
    "* Describe and fit memoryless linear models\n",
    "* Describe and fit linear time-series models with delays\n",
    "* Fit linear models with multiple target outputs\n",
    "* Select the optimal delay via cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the data\n",
    "\n",
    "The data in this lab comes from neural recordings described in:\n",
    "\n",
    "<a href=\"http://jn.physiology.org/content/106/2/764.short\">\n",
    "Stevenson, Ian H., et al. \"Statistical assessment of the stability of neural movement representations.\" Journal of neurophysiology 106.2 (2011): 764-774</a>\n",
    "\n",
    "Neurons are the basic information processing units in the brain.  Neurons communicate with one another via *spikes* or *action potentials* which are brief events where voltage in the neuron rapidly rises then falls.  These spikes trigger the electro-chemical signals between one neuron and another.  In this experiment, the spikes were recorded from 196 neurons in the primary motor cortex (M1) of a monkey using an electrode array implanted onto the surface of a monkey's brain.  During the recording, the monkey performed several reaching tasks and the position and velocity of the hand was recorded as well.  \n",
    "\n",
    "The goal of the experiment is to try to *read the monkey's brain*:  That is, predict the hand motion from the neural signals from the motor cortex.\n",
    "\n",
    "We first load the key packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import r2_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full data is available on the CRCNS website  http://crcns.org/data-sets/movements/dream.  This website has a large number of great datasets and can be used for projects as well.  However, the raw data files can be quite large.  To make the lab easier, the [Kording lab](http://kordinglab.com/) at UPenn has put together an excellent [repository](https://github.com/KordingLab/Neural_Decoding) where they have created simple pre-processed versions of the data.  You can download the file `example_data_s1.pickle` from the [Dropbox link](https://www.dropbox.com/sh/n4924ipcfjqc0t6/AADOv9JYMUBK1tlg9P71gSSra/example_data_s1.pickle?dl=0).  Alternatively, you can directly run the following code.  This may take a little while to download since the file is 26 MB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fn_src = 'https://www.dropbox.com/sh/n4924ipcfjqc0t6/AADOv9JYMUBK1tlg9P71gSSra/example_data_s1.pickle?dl=1'\n",
    "fn_dst = 'example_data_s1.pickle'\n",
    "\n",
    "import os\n",
    "from six.moves import urllib\n",
    "\n",
    "if os.path.isfile(fn_dst):\n",
    "    print('File %s is already downloaded' % fn_dst)\n",
    "else:        \n",
    "    urllib.request.urlretrieve(fn_src, fn_dst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file is a *pickle* data structure, which is a package to serialize python objects into data files.  Once you have downloaded the file, you can run the following command to retrieve the data from the pickle file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('example_data_s1.pickle', 'rb') as fp:\n",
    "    X,y = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix `X` is matrix of spike counts where `X[i,j]` is the number of spikes from neuron `j` in time bin `i`.\n",
    "The matrix `y` has two columns:\n",
    "* `y[i,0] = ` velocity of the monkey's hand in the x-direction\n",
    "* `y[i,1] = ` velocity of the monkey's hand in the y-direction\n",
    "Our goal will be to predict `y` from `X`.  \n",
    "\n",
    "Each time bin represent `tsamp=0.05` seconds of time.  Using `X.shape` and `y.shape` compute and print:\n",
    "* `nt = ` the total number of time bins\n",
    "* `nneuron = ` the total number of neurons\n",
    "* `nout = ` the total number of output variables to track = number of columns in `y`\n",
    "* `ttotal = ` total time of the experiment is seconds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsamp = 0.05  # sampling time in seconds\n",
    "\n",
    "# TODO\n",
    "# nt = ...\n",
    "# nneuron = ..\n",
    "# nout = ...\n",
    "# ttotal = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a Memoryless Linear Model\n",
    "\n",
    "Let's first try a simple linear regression model to fit the data.\n",
    "\n",
    "First, use the `train_test_split` function to split the data into training and test.  Let `Xtr,ytr` be the training data set and `Xts,yts` be the test data set.  Use `test_size=0.33` so `1/3` of the data is used for test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# TODO\n",
    "# Xtr, Xts, ytr, yts = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, fit a linear model using `Xtr,ytr`.  Make a prediction  `yhat` using `Xts`.  Compare `yhat` to `yts` to measure `rsq`, the  `R^2`.  You can use the `r2_score` method. Print the `rsq` value.  You should get `rsq` of around `0.45`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# yhat = ...\n",
    "# rsq = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is useful to plot the predicted vs. true values.  Since we have two outputs, create two `subplots` using the `plt.subplot()` command.  In plot `i=0,1`, plot `yhat[:,i]` vs. `yts[:,i]` with a scatter plot.   Label the axes of the plots.  You may also use the command:\n",
    "\n",
    "    plt.figure(figsize=(10,5))\n",
    "    \n",
    "to make the figures a little larger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting Models with Delay\n",
    "\n",
    "One way we can improve the model accuracy is to used delayed version of the features.  Specifically, the model we used above mapped the features \n",
    "\n",
    "    yhat[i,k] = \\sum_{j=0}^{p-1} X[i,j]*w[j,k] + b[k]\n",
    "    \n",
    "where `p` is the number of features and `w[j,k]` is a matrix of coefficients.  In this model,  `yhat[i,:]` at time `i` was only dependent on the inputs  `X[i,:]` at time `i`.  In signal processing, this is called a *memoryless* model.  However, in many physical systems, such as those that arise in neuroscience, there is a delay between the inputs `X[i,:]` and the outputs `y[i]`.  For such cases, we can use a model of the form,\n",
    "\n",
    "    yhat[i+d,k] = \\sum_{j=0}^{p-1} \\sum_{m=0}^d X[i+m,j]*W[j,m,k] + b[k]\n",
    "    \n",
    "where `W` is a 3-dim array of coefficients where:\n",
    "\n",
    "    W[j,m,k] is the influence of the input X[i+m,j] onto output y[i+d,k]\n",
    "\n",
    "\n",
    "In signal processing, this model is called an *FIR* filter and `W[j,:,k]` is the *impulse response* from the `j`-th input to the `k`-th output.  The point is that the output at time `i+d` depends on the inputs at times `i,i+1,...,i+d`.  Hence, it depends on the last `d+1` time steps, not just the most recent time.\n",
    "\n",
    "To translate this into a linear regression problem, complete the following function that creates a new feature and target matrix where:\n",
    "\n",
    "    Xdly[i,:] has the rows X[i,:], X[i++1,:], ..., X[i+dly,:]\n",
    "    ydly[i,:] = y[i+dly,:]\n",
    "    \n",
    "Thus, `Xdly[i,:]` contains all the delayed fetaures for the target `yhat`.  Note that if `X` is `n x p` then `Xdly` will be `n-dly x (dly+1)*p`.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_dly_data(X,y,dly):\n",
    "    \"\"\"\n",
    "    Create delayed data\n",
    "    \"\"\"    \n",
    "    # TODO\n",
    "    # Xdly = ...\n",
    "    # ydly = ...\n",
    "    \n",
    "    return Xdly, ydly\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now fit an linear delayed model with `dly=6` additional delay lags.  That is,\n",
    "* Create delayed data `Xdly,ydly=create_dly_data(X,y,dly=6)`\n",
    "* Split the data into training and test as before\n",
    "* Fit the model on the training data\n",
    "* Measure the `R^2` score on the test data\n",
    "\n",
    "If you did this correctly, you should get a new `R^2` score around 0.69.  This is significantly better than the memoryless models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the predicted vs. true values as before. You should visually see a better fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note*:  Fitting an FIR model with the above method is very inefficient when the number of delays, `dly` is large.  In the above method, the number of columns of `X` grows from `p` to `(dly+1)*p` and the computations become expensive with `dly` is large.  We will describe a much faster way to fit such models using gradient descent when we talk about convolutional neural networks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting the Optimal Delay via Model Order Selection\n",
    "\n",
    "In the previous example, we fixed `dly=6`.  We can now select the optimal delay using model order selection.  Since we have a large number of data samples, it turns out that the optimal model order uses a very high delay.  Using the above fitting method, the computations take too long.  So, to simplify the lab, we will first just pretent that we have a very limited data set.\n",
    "\n",
    "Compute `Xred` and `yred` by taking the first `nred=6000` samples of the data `X` and `y`.  This is about 10% of the overall data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nred = 6000\n",
    "\n",
    "# TODO\n",
    "# Xred = ...\n",
    "# yred = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now complete the following code to implement K-fold cross validation with `nfold=5` and values of delays `dtest = [0,1,...,dmax]`.  \n",
    "\n",
    "The code also includes a progress bar using the `tqdm` package.  This is very useful when you have a long computation.\n",
    "\n",
    "Note:  Some students appeared to use the `mse` metric (i.e. RSS per sample) instead of `R^2`.  That is fine.  For the solution, I have computed both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import  sklearn.model_selection \n",
    "import tqdm.notebook\n",
    "\n",
    "nfold = 5  # Number of folds\n",
    "dmax = 15  # maximum number of delays\n",
    "\n",
    "# TODO:  Create a k-fold object\n",
    "# kf = sklearn.model_selection.KFold(...)\n",
    "\n",
    "# TODO:  Model orders to be tested\n",
    "# dtest = vector [0,1,...,dmax]  \n",
    "# nd = len(dtest)\n",
    "\n",
    "# TODO.  \n",
    "# Initialize a matrix Rsq to hold values of the R^2 across the model orders and folds.\n",
    "# Alternatively, you can also create an RSS matrix\n",
    "\n",
    "# Create a progress bar.  Note there are nd*nfold total fits.\n",
    "pbar = tqdm.notebook.tqdm(\n",
    "    total=nfold*nd, initial=0,\n",
    "    unit='fits', unit_divisor=nd, desc='Model order test')\n",
    "\n",
    "\n",
    "for it, d in enumerate(dtest):\n",
    "    # TODO:  \n",
    "    # Create the delayed data using the create_dly_function from the reduced\n",
    "    # data Xred, yred\n",
    "    #     Xdly, ydly = create_dly_data(...)    \n",
    "        \n",
    "    # Loop over the folds\n",
    "    for isplit, Ind in enumerate(kf.split(Xdly)):\n",
    "\n",
    "        # Get the training data in the split\n",
    "        Itr, Its = Ind        \n",
    "\n",
    "        # TODO\n",
    "        # Split the data (Xdly,ydly) into training and test\n",
    "        # Xtr = ...\n",
    "        # ytr = ...\n",
    "        # Xts = ...\n",
    "        # yts = ...\n",
    "        \n",
    "        # TODO:  Fit data on training data        \n",
    "\n",
    "        # TODO:  Measure the R^2 vale on test data and store in the matrix Rsq\n",
    "        \n",
    "        pbar.update(1)\n",
    "pbar.close()       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the mean and standard error of the `R^2` values as a function of the model order `d`.  Use a `plt.errorbar` plot.  Label your axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the optimal order `d` with the normal rule (i.e. highest test `R^2`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now find the optimal model order via the one SE rule (i.e. highest test `R^2` within on SE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}