05-Neuroanalysis with Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neuroanalysis with Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although you may not believe it, you have now learned enough NumPy and Matplotlib to **tackle almost any problem in neuroanalysis!!** But what exactly do I mean with neuroanalysis? In this session and the next, we'll take a look at some of the common kinds of problems and anslysis you will encounter when working with neural data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First thing first, let's go ahead and import essential packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with neuronal recordings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you perform electrophysiology recordings on a neuron or perhaps a population of neurons, you will be working with a recording of a **neuron's electrical potential**. This might be made available to you in varioud **data formats**. When you want to work with them, usually the very first question is: **how do I load the data?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading data files into Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, when you worked with NumPy arrays, most of the data was generated by you, the programmer. With use of functions like `input` we also saw how you could prompt the **user** of your code to enter some data into your code, usually by entering data one at a time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, in scientific computing, you would very likely want to work with **data that has been generated elsewhere**. The source of data can be an eletrophysiology system, microscopy imaging system, eye tracker, or even surveys filled out by human subjects. In modern day computing, these data are then typically stored and transported as **files** with some data format. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, if you want to get data, you will have to somehow load data from files into Python. There are many ways to do this, and we will take a look at few options you have."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading simple textual files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the simplest file format is **plain-text** files. As the name suggests, these are files that only contain textual characters inside. For an example, let's take a look at `recordings.dat` file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are going to use Jupyter's `%load` magic to open up the content of the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load data/mouse.dat\n",
    "m122\n",
    "m104\n",
    "m48\n",
    "m690\n",
    "m329"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to **load the content of the file**. To do so, we are going to use the `open` function in Python to `open` a file!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/mouse.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This returns a **file handle** - a special kind of object that **represents a file**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the entire content\n",
    "content = f.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run read for the second time\n",
    "content_again = f.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "content_again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you are done using a file, you should always **close** the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.close() # closes the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.read() # operation on a closed file triggers an error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### On a tangent: Escape sequences in a string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes you would want to express **special key** results such as the `Enter` key resulting in a **line break** and `Tab` key resulting in an indentation. In Python (and in fact almost every other languages), you express these with **escape sequences**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a string with line break `\\n`\n",
    "message = 'Hi! This message\\nshould appear over\\nthree lines!'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Escape sequences like `\\n` will not be treated too specially if you just look into the content of a string. The effect really kicks in when you **print** the string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(message)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('This\\tshould\\tappear\\ttab\\tshifted\\nOccuring\\tover\\ttwo\\tlines')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interesting escape sequence is **a single character** as far as Python is concerned:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var1 = 'abc'\n",
    "print(var1)\n",
    "print(len(var1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var2 = 'a\\nc'\n",
    "print(var2)\n",
    "print(len(var1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are [many other escape sequences](https://docs.python.org/3/reference/lexical_analysis.html#string-and-bytes-literals), the newline character `\\n` and tab character `\\t` are by far the two most common escape sequences you would encounter and thus would want to know!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Escaping characters with special meaning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Escape sequence can be useful if you want to use `'` and `\"` inside a string literal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While you could use `'` and `\"` inside a string literal as long as you use the other quote kind:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('This method \"works\" well')\n",
    "print(\"But limits your 'choice'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Now I can put \"double quotes\" and \\'single quotes\\' in a same string!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading file content line by line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running `read` method on the file handle gives you the entire file content in one string. Although this is a start, this may not be the most convenient way to work with a file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thankfully, file handle is an **iterable** object - meaning you can iterate through the file content!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/mouse.dat')\n",
    "\n",
    "for line in f:\n",
    "    print(line)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that there is an extra line between each because there is `\\n` character at the end of each line that is read."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can remove any leading and trailing spaces and new line characters -- collectively referred to as **white space characters** -- by using `strip` method on the string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/mouse.dat')\n",
    "\n",
    "for line in f:\n",
    "    print(line.strip())  # strips leading/trailing whiltespace characters\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can of course store the line content into a list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/mouse.dat')\n",
    "mouse_list = []\n",
    "for line in f:\n",
    "    mouse_list.append(line.strip())  # strips leading/trailing whiltespace characters\n",
    "    \n",
    "f.close()\n",
    "    \n",
    "mouse_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, it turns out that file object has a method called `readlines` that **almost** does what we want!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.readlines?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/mouse.dat')\n",
    "mouse_list = f.readlines()\n",
    "\n",
    "mouse_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only issue is that each data entry has the trailing `\\n` on it. Although you can do a full fledged for loop, this is where **list comprehension** can be very handy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_mouse_list = [m.strip() for m in mouse_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_mouse_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading comma separated values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another common file format you might encounter is a textual file where **distinct values are separated by comma**. These files are typically referred to as **comma-separated values** and have file extension **csv**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load data/position.csv\n",
    "0.65589,0.57823,0.33661,0.16374,0.12180,-0.08730,-0.32791,-0.02155,-0.56670,-0.70968,-0.83033,-0.94229,-0.62620,-0.83811,-0.84556,-1.00000,-0.99296,-0.99012,-0.87674,-0.86348,-0.92978,-0.93794,-0.74789,-0.54437,-0.69714,-0.57095,-0.47770,-0.05630,-0.05112,0.02798,0.15682,0.56518,0.42280,0.54062,0.73925,0.46023,0.93578,0.83758,0.92696,0.83959,0.98542,0.93230,0.93887,0.98559,0.91488,0.99161,0.66137,0.93301,0.55899,0.32709,0.71629,0.38173,0.28290,0.27180,-0.27958,-0.36215,-0.28922,-0.49119,-0.68146,-0.32334,-0.60336,-0.90789,-0.89853,-0.99541,-0.96631,-0.99483,-0.99976,-0.99998,-0.95507,-0.84919,-0.83438,-0.91877,-0.85928,-0.65672,-0.54651,-0.53178,-0.40231,-0.01891,-0.15738,-0.10022,0.26732,0.25906,0.66773,0.56343,0.61935,0.80379,0.87985,0.75094,0.84290,0.98704,0.97718,0.98602,0.88545,0.99304,0.84986,0.88180,0.74419,0.71700,0.58715,0.54574"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/position.csv')\n",
    "data = f.read()\n",
    "f.close()\n",
    "\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So problem is that now what you have is a big string of values separated by comma. You would first want to **split them at comma**. Well, that's what the `split` method is for!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "split_data = data.split(',')  # split at comma\n",
    "\n",
    "split_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the data is still in strings, so we want to convert each of these values into `float`. Let's use list comprehension again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# converts every element of split_data into float\n",
    "float_data = [float(v) for v in split_data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "float_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quickly plot this data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(float_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like the subject was moving back and forth!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### On the tangent: plotting data with `plot`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have been surprised to see that Matplotlib's `plot` accepts a plain list instead of a NumPy array and that it was okay with just a single list. Somehow, it didn't need x values!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plot` actually happily accepts any **iterable** data as its x or y value. So you can use `plot` with plain lists!\n",
    "\n",
    "Also, if you don't supply two lists (for x and y) and just give one list, then the list is automatically considered to be the y-values. `plot` will silently prepare it's own x values that counts up the data points from 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Letting NumPy do the hard work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, you may be surprised how much work is involved in just loading such simple files! Indeed, **figuring out how to load the data could be the hardest part of your project** depending on the project and how data is stored and passed around!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortunately, packages like NumPy comes with functions to make loading some files less painful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading CSV files with NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy comes with `loadtxt` function that can work with a large variety of file formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.loadtxt?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.loadtxt('data/position.csv', delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It *could* deal with more complex file content."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.loadtxt('data/block_data.csv', skiprows=1, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, NumPy `loadtxt` starts to show limitation when the multiple **data types** are found in a single file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.loadtxt('data/named_data.csv', skiprows=1, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When loading CSV or even excel spread sheet data, you are better off using more **table oriented packages** such as `pandas`. Just a quick demo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/named_data.csv')\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[1:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with NumPy specific file format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although you could use formats like **line separated values** and **commad separated values** to store 1-D and 2-D data, saving and loading higher dimensional data can be a challenge. Fortunately, NumPy comes with its own data format that you can use to store and load NumPy arrays conveniently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving in NumPy file format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a 4d array filled with random numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.random.randn(10, 5, 3, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can **save this data into a NumPy data file** using `np.save` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('data/random4d', data) # by default puts `.npy` extension"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it! You have just saved your random array into a file called `data/random4d.npy`!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading in NumPy file format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, saving is just half of the story - we must be able to load it! As you might have guessed, you would use `np.load` function to load data saved in NPY file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loaded_data = np.load('data/random4d.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loaded_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use `np.all` to verify that our `loaded_data` is the same as the saved `data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(data == loaded_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But notice that these two are **separate copies**: i.e. changing one does **not** change the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:] = 0 # fill the entire data with 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loaded_data # loaded copy is unaffected!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have seen various ways of loading files and even learned to save your NumPy array into file, we are ready to work with neural data! But before we move on, let's visit the topic of saving to a file quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Before moving on: saving to a file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might have noticed that I haven't really covered how to **save data into a textual file**! This is because:\n",
    "\n",
    "1. You are far more likely to want to load textual files then to write one out from Python\n",
    "2. More often than not, there are better tools available for saving your data out to one of many common data formats rather than writing textual files by yourself!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, I'm making a deliberate choice of **not covering how to write your own textual file** here. For interested individuals, I have prepared a separate notebook that covers not only writing textual files, but also on reading and writing binary files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading neuronal data provided to you"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have learned quite a bit about files, let's go ahead and load **a sample** neuronal trace provided to you as a NumPy file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity = np.load('data/neural_trace.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity # careful! this will printout a long list!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the loaded NumPy array is a rather simple 1D array with a lot of data point!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now spend some time analyzing this data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing your data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the very first thing you should do when you get a dataset is to **visualize** them! Doing so will often gives you much deeper understanding of your data, and also often given you a chance to detect any anomaly in your data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you are told that the data represents **a time series of a neuron's electrical activity**. So let's go ahead and plot it as a time series with `plot`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(activity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so this does indeed look like an electricial activity of a **spiking neuron**, with sharp rises corresponding to some spikes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the data looks quite **noisy** so that every spike is not a picture perfect action potential you may have come to expect from textbooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a closer look at a spike by **slicing** the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(activity[45000:46000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we actually don't know what x-axis represents. It's is very likely related to time, but we don't know what each point represents.\n",
    "\n",
    "Say now you are told that each point represents a value of the neuron's electric activity read out at **32kHz** - that is, 32,000 samples are taken every second! You are now also told that the activity is measured in mV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's prepare a suitable x values. We can convert the sample number into a time value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_number = np.arange(len(activity)) # count up for each sample\n",
    "sample_time = sample_number / 32000 # convert into time in seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(sample_time, activity)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Membrane potential (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing statistics on the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have visualized and have some idea about the data, it's time to start analyzing them! Often, it is a good idea to compute some **statistics** on your data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean, median, and standard deviation of the electric activity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go ahead and compute mean, median, and standard deviation of the electric activity. All of these are available as either NumPy function or methods on the NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(activity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(activity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.median(activity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, it can be useful to look at the **distribution** of the data itself. A histogram can give you a sense of variability in your dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(dpi=300)\n",
    "\n",
    "plt.hist(activity)\n",
    "\n",
    "plt.xlabel('Membrane potential (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Default setting on the `hist` plot gives rather crude histogram. You can use **much finer binning** by specifying the number of bins to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(dpi=300)\n",
    "\n",
    "plt.hist(activity, bins=100) # use more bins to describe the distribution\n",
    "\n",
    "plt.xlabel('Membrane potential (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if you look at the data, there seems to be two *segments* to the activity - in the first half there is not much activity but there is a lot more in the latter half."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go ahead and split the data into half and analyze them separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the index of the half point\n",
    "half_point = len(activity) // 2\n",
    "\n",
    "half_point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_half = activity[:half_point]\n",
    "second_half = activity[half_point:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(len(first_half)) / 32000\n",
    "plt.plot(t, first_half)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(len(second_half)) / 32000\n",
    "plt.plot(t, second_half)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute the statistics! Let's write a utility function to compute and print out stats for the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_stats(data):\n",
    "    print('Mean:', data.mean())\n",
    "    print('Std:', data.std())\n",
    "    print('Median:', np.median(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('For first half')\n",
    "print_stats(first_half)\n",
    "print('---------------')\n",
    "print('For second half')\n",
    "print_stats(second_half)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the results, we can see that the standard deviation is lower for the first half. We can consider the first half to be the **resting activity** of the neuron, and the standard deviation during this period is often considered the **baseline noise standard deviation**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Characterizing neural responses to stimulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of working with so-called raw electric recordings, let's assume that you work with some measure of neural activity such as **firing rates** or **total spike counts** for each **trial** in an experiment where you present to the animal some **stimulus**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say that we have performed an experiment where you are recording from a single neuron in the primary visual cortex (V1) of mouse, and on each trial, you present to the mouse a visual stimulus (e.g. an image of a racoon) at **different contrast** to the mouse:\n",
    "\n",
    "![contrast](images/contrasts.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the experiment, we got a NPY file containing the data about the contrast of the image presented, and the corresponding acitivity of the neuron measured in firing rate in a file called `data/contrast_exp.npy`. Let's go ahead and load it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.load('data/contrast_exp.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So it is a 2D array of shape 2 x 55. Let's assume we know that the two rows correspond to the contrast and the neuron's response in firing rates.\n",
    "\n",
    "Again, the very first thing you should try is **visualization**!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = data[0]\n",
    "y = data[1]\n",
    "\n",
    "fig = plt.figure(dpi=300)\n",
    "plt.plot(x, y)\n",
    "\n",
    "plt.xlabel('Contrast')\n",
    "plt.ylabel('Firing rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like a mess because this is a data where data points are organized in any specific sequence! For these kinds of data, `scatter` plot would be a better fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrasts = data[0]\n",
    "firing_rates = data[1]\n",
    "\n",
    "fig = plt.figure(dpi=300)\n",
    "plt.scatter(contrasts, firing_rates)\n",
    "\n",
    "plt.xlabel('Contrast')\n",
    "plt.ylabel('Firing rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better! Just by looking at the plot, we can tell that the firing rate of the neuron appears to be increase with the contrast of the image. In fact, it almost looks like it follows a straight line!\n",
    "\n",
    "Whenever you sense that there is a **linear relationship** between two variables - here contrast and neuron's firing rate, it is a good idea to perform **linear regression** on your dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear regression actually falls into a category of analysis called **model fitting** where you try to explain relationships in your data with **a model**. Specfically, in linear regression, you are assuming that **there is some linear relationship** between the two variables, in this case, the contrast and the firing rate. \n",
    "\n",
    "In other words, you are assuming that the **underlying model is**:\n",
    "\n",
    "$ \\text{firing rate} = \\text{slope} * \\text{contrast} + \\text{intercept}$\n",
    "\n",
    "For some value of slope and intercept!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if you think about it, this is not possible, because the contrast and firing rate **doesn't seem to be following a perfect line**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get around this issue by assuming that **your data is noisy**! That is, for every measurement you make in firing rate, there is some noise that you cannot control!\n",
    "\n",
    "$ \\text{firing rate} = \\text{slope} * \\text{contrast} + \\text{intercept} + \\text{noise}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under the standard linear regression, you assume that your noise has a **normal distribution** that is, if you were to plot the histogram of noise, it should look something like:\n",
    "\n",
    "![noise distribution](images/normal.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I will show you two ways of performing linear regression. First through the use of `scipy` package's function `linregress`. Second, we will look at how to perform linear regression entirely ourselves mathematically!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SciPy package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus far we have talked about two major scientific packages in Python: NumPy and Matplotlib. Now I'm introducing you another one called SciPy, standing for Scientific Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike NumPy and Matplotlib that had more or less specific purpose (numerical array and plotting, respectively), SciPy is really a large collection of various packages each geared towards different aspects of scientific computations. Here we are going to take a look at one of its **subpackage** called `stats`, which is about **statistical analysis and tests**. Because SciPy is such a big package, you typically do not import the package itself, but rather import specific subpackage and its content that you want to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get a sense of what comes with SciPy by visiting [their website](https://docs.scipy.org/doc/scipy/reference/tutorial/general.html)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here let's go ahead and import the `stats` subpackage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go ahead and take some time to admire the number of functions and properties inside this package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's take a look at the function we are ultimately interested in: `linregress` and look at its documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.linregress?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It gives quite descriptive documentation complete with a usage example. Let's now use this function to fit a linear model to our data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = stats.linregress(contrasts, firing_rates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it! Let's now take a look at our result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_value, p_value, std_err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have values for the slope and intercept, let's use our model to compute the **expected value of firing rate** for each contrast, and plot them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some values of contrasts to make prediction for\n",
    "x = np.linspace(0, 100, 100)\n",
    "\n",
    "# predicted firing rates\n",
    "y = x * slope + intercept\n",
    "\n",
    "plt.plot(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In itself, this is not too interesting - obviously we want to compare this model **against** the real data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some values of contrasts to make prediction for\n",
    "x = np.linspace(0, 100, 100)\n",
    "\n",
    "# predicted firing rates\n",
    "y = x * slope + intercept\n",
    "\n",
    "plt.plot(x, y)\n",
    "\n",
    "# plot the actual data on top\n",
    "plt.scatter(contrasts, firing_rates, color='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks like a reasonable fit! Not, let's actually compute the expected firing rate for each value of contrast found in the experiment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_firing_rates = contrasts * slope + intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and compute the **difference** between the expected values and the actually observed values of the firing rates!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "difference = firing_rates - expected_firing_rates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we want to take a look at the **distribution of the deviations, or noise**!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(difference, bins=10)\n",
    "plt.xlabel('Difference in firing rates')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although not perfect, the difference appears to be reasonably normally distributed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Doing linear regression on your own"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While you can use implementation of linear regression as provided by packages like SciPy, you can actually compute and perform linear regression on your own! This does involve some math involving vectors and matrices. I'm including the steps strictly for the interested individuals:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For linear regression with intercept, this is the same as solving the expression:\n",
    "\n",
    "$$\n",
    "Y = Xw\n",
    "$$\n",
    "\n",
    "where $X$ is an $n \\times 2$ matrix containing all x data as first column and a vector of 1s in the second column, and $Y$ is $n \\times 1$ vector containing y data. $w$ is a element 2 vector whose first element is the slope and the second is the intercept of linear regression. In such case, the solution can be shown to be:\n",
    "\n",
    "$$\n",
    "w = (X^\\top X)^{-1}X^\\top Y\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get started in computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = contrasts # for easier manipulation\n",
    "y = firing_rates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first prepare $X$. Recall that it's first column is the x data and the second column is filled with 1. To get this, we construct a vector of ones and **stack** the two vectors into a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ones = np.ones(x.shape) # prepare a vector of ones with same length as x\n",
    "\n",
    "X = np.stack([x, ones], axis=1) # stack x and vector of ones to get X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we prepare $Y$ by adding an extra dimension to `y`, making it a column vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y =  y[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we get to compute $w$. The equation involves **inverting a matrix**. This is a special operation that needs a function which is provided by NumPy under it's `linalg` subpackage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.linalg import inv # import matrix inversion function used to take inverse of a matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now go ahead and compute $w$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = inv(X.T @ X) @ X.T @ Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at `w`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same exact value provided by `scipy.stats.linregress` function!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below I step into more advanced analysis for those of you interested."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing data for further analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have visualized the data and computed some statistics, we have pretty good sense of how the data looks. Often times, you would want to **process** your data further before you perform additional analysis. For example, when working with neuronal electricial activity, you'd be interested in **detecting spikes** and extract **when the spike occurred**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, let's load the data again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity = np.load('data/neural_trace.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spike detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Developing spike detection algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we set out, let me state that **spike detection** is a very in-depth topic, with a lot of research involved in how to **extract spikes** from neuronal recordings of electric activity. Given this, what we are about to implement is just about the simplest spike detection algorithm out there.\n",
    "\n",
    "However, it serves very well for the purpose of illustration on how in principle spike detection would work, and more importantly, how to use Python in achieving non-trivial computations!\n",
    "\n",
    "With that out of the way, let's continue!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic idea behind spike detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look again at a single spike trace in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_number = np.arange(len(activity)) # count up for each sample\n",
    "t = sample_number / 32000 # convert into time in seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "segment = np.logical_and(t > 1.36, t < 1.38)\n",
    "\n",
    "plt.plot(t[segment], activity[segment])\n",
    "\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Membrane potential (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a spike (an action potential) occurs, the membrane potential (the electric activity) quickly rises. Perhaps we can try to detect when the potential goes above certain **threshold**. To develop an algorithm, let's go ahead and focus on this particular segment of activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = t[segment]\n",
    "potential = activity[segment]\n",
    "\n",
    "plt.plot(ts, potential)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try to detect when the potential is above a certain value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = -40\n",
    "\n",
    "above_thr = potential > threshold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives back an array of `True`'s and `False`'s. It turns out that, you can convert a boolean array into float, and when you do so, it will convert `True`'s into `1`'s and `False`'s into `0`'s. You can convert types of array with `.astype` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ts, above_thr.astype(float))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, if you straight up plot a boolean array, `plot` will do this conversion for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ts, above_thr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, it would be nice to see this side by side with the membrane potential. Let's try plotting them on the same axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ts, potential)\n",
    "plt.plot(ts, above_thr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting them on the same axes doesn't quite give good visibility due to scaling difference. Let's use *subplots*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = -40\n",
    "\n",
    "above_thr = potential > threshold\n",
    "\n",
    "fig = plt.figure(dpi=300)\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(ts, potential)\n",
    "plt.axhline(threshold, color='orange')\n",
    "\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(ts, above_thr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, that looks like a good start. We now would like to know **when did the potential go above the threshold**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with toy example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's think on how to do this with a much simpler data. Some like below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = np.array([0, 0, 0, 1, 1, 1, 0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, we know that it **goes above threshold on index 3** and comes back down on **index 7**. The question is, how can we detect this without manually looking at it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You could imagine writing a for-loop that finds when the value went from 0 to 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_value = 0\n",
    "for i, v in enumerate(sample):\n",
    "    if last_value == 0 and v == 1:\n",
    "        print('Changed at index {}!'.format(i))\n",
    "    last_value = v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it is often better to try to avoid for-loops: if something can be done all withint NumPy, it is often significantly faster.\n",
    "\n",
    "Here, it would be nice if we can figure out the change in value from one index to another. Well, it turns out that you can use `np.diff`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.diff(sample) # gives adjacent value difference\n",
    "\n",
    "delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice! Now we can see that whereever it changed from 0 to 1 has difference of 1. In fact, we can also detect whereever it changed from 1 to 0 as difference is -1!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where(delta == 1) # find out where difference was 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Funnily enough, we are told that the index is at 2, not like 3 that we expected. This is because `np.diff` caused the result to shrink by 1 in size! To get around this, we can **pad** `delta` with a 0 at the beginning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "padded_delta = np.pad(delta, (1, 0), mode='constant', constant_values=0)\n",
    "\n",
    "padded_delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the value 1 occurs at the exact index where value went from 0 to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where(padded_delta == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying back to real data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we know how to solve the problem of **when did the potential go above the threshold** on a toy example, let's apply it back to the real problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we had `above_thr` that contained `True`'s and `False`'s for where the potential was above the threshold. Let's use the fact that you can treat the array of boolean as array of 0s and 1s, and just apply the algorithm we developed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.diff(above_thr.astype(float)) # find out changes in values\n",
    "padded_delta = np.pad(delta, (1, 0), mode='constant', constant_values=0)\n",
    "\n",
    "np.where(padded_delta == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = -40\n",
    "\n",
    "above_thr = potential > threshold\n",
    "\n",
    "delta = np.diff(above_thr.astype(float)) # find out changes in values\n",
    "padded_delta = np.pad(delta, (1, 0), mode='constant', constant_values=0)\n",
    "\n",
    "spike_pos = padded_delta == 1\n",
    "spike_time = ts[spike_pos]\n",
    "\n",
    "fig = plt.figure(dpi=300)\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(ts, potential)\n",
    "plt.axhline(threshold, color='orange')\n",
    "# draw vertical line at detected spike position\n",
    "plt.axvline(spike_time, linestyle='--', color='red')\n",
    "\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(ts, above_thr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying the algorithm to the entire dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have developed the algorithm on a small dataset, it's time to apply it to the original dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = -40\n",
    "\n",
    "above_thr = activity > threshold\n",
    "\n",
    "delta = np.diff(above_thr.astype(float)) # find out changes in values\n",
    "padded_delta = np.pad(delta, (1, 0), mode='constant', constant_values=0)\n",
    "\n",
    "spike_pos = padded_delta == 1\n",
    "spike_times = t[spike_pos] # this will be multiple values\n",
    "\n",
    "fig = plt.figure(dpi=300)\n",
    "\n",
    "plt.subplot(3, 1, 1)\n",
    "plt.plot(t, activity)\n",
    "plt.axhline(threshold, color='orange')\n",
    "\n",
    "plt.subplot(3, 1, 2)\n",
    "plt.plot(t, above_thr)\n",
    "\n",
    "\n",
    "plt.subplot(3, 1, 3)\n",
    "plt.plot(t, activity)\n",
    "for spike_time in spike_times:\n",
    "    plt.axvline(spike_time, linestyle='--', color='red')"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}