{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 04 - Handling spiking data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals\n",
    "\n",
    "- Describe spike trains through binning and spike density functions\n",
    "- Describe spike trains by their interspike intervals (ISIs)\n",
    "- Compute an autocorrelation function (ACF) and crosscorrelation function (CCF)\n",
    "- Generating fake spike data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine firing rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:10:05.101891Z",
     "start_time": "2017-07-26T18:10:04.148685Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "%matplotlib inline\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import nept\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal\n",
    "\n",
    "# define where your data folder is located\n",
    "data_path = os.path.join(os.path.abspath('.'), 'data')\n",
    "data_folder = os.path.join(data_path, 'R042-2013-08-18')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:10:13.564197Z",
     "start_time": "2017-07-26T18:10:13.555209Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load the info file, which contains experiment-specific information\n",
    "sys.path.append(data_folder)\n",
    "import r042d3 as info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:10:15.172102Z",
     "start_time": "2017-07-26T18:10:15.138596Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load spikes (.t and ._t) from this experiment\n",
    "spikes = nept.load_spikes(data_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:12:40.105124Z",
     "start_time": "2017-07-26T18:12:40.098620Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's limit our investigation to one neuron\n",
    "neuron_idx = 31\n",
    "these_spikes = spikes[neuron_idx]\n",
    "\n",
    "# And restrict the time so we can more easily see what's going on\n",
    "start = 2500.0\n",
    "stop = 2700.0\n",
    "filtered_spikes = these_spikes.time_slice(start, stop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:12:40.896727Z",
     "start_time": "2017-07-26T18:12:40.747620Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the spikes\n",
    "plt.plot(filtered_spikes.time, np.ones(len(filtered_spikes.time)), '|', ms=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:12:30.834267Z",
     "start_time": "2017-07-26T18:12:30.683161Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(these_spikes.time, np.ones(len(these_spikes.time)), '|', ms=30)\n",
    "plt.xlim(2500, 2700)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:13:11.127231Z",
     "start_time": "2017-07-26T18:13:11.122726Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create an AnalogSignal used to define the time edges for the binned spikes\n",
    "edges = nept.AnalogSignal(np.ones(20), np.linspace(start, stop, 20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:17:58.407148Z",
     "start_time": "2017-07-26T18:17:58.402645Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Bin the spikes\n",
    "window_advance = 0.5\n",
    "time_edges = nept.get_edges(edges, window_advance, lastbin=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:21:41.745201Z",
     "start_time": "2017-07-26T18:21:41.579559Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the spikes\n",
    "plt.plot(filtered_spikes.time, np.ones(len(filtered_spikes.time))-1.5, '|', ms=10)\n",
    "\n",
    "# Plot the number of spikes in each bin\n",
    "plt.hist(filtered_spikes.time, time_edges, histtype='step')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:17:31.744260Z",
     "start_time": "2017-07-26T18:17:31.691192Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now let's look at a spike density function (SDF; convolved spike train)\n",
    "# for this SDF we need a smaller bin size\n",
    "bin_size = 0.001\n",
    "\n",
    "sdf_edges = np.arange(start, stop, bin_size)\n",
    "sdf_centers = sdf_edges[:-1]+bin_size/2\n",
    "\n",
    "# Make a gaussian filter\n",
    "gaussian_window = 1.0 / bin_size\n",
    "gaussian_std = 0.02 / bin_size\n",
    "\n",
    "gaussian_kernel = scipy.signal.gaussian(gaussian_window, gaussian_std)\n",
    "gaussian_kernel /= np.sum(gaussian_kernel)\n",
    "gaussian_kernel /= bin_size\n",
    "\n",
    "# Bin the spikes\n",
    "spike_count = np.histogram(filtered_spikes.time, bins=sdf_edges)[0]\n",
    "\n",
    "# Convolve the binned spikes by the gaussian filter\n",
    "convolved_spiketimes = scipy.signal.convolve(spike_count, gaussian_kernel, mode='same')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:21:26.504260Z",
     "start_time": "2017-07-26T18:21:26.307135Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the spikes\n",
    "plt.plot(filtered_spikes.time, np.ones(len(filtered_spikes.time))-5, '|', ms=10)\n",
    "\n",
    "# Plot the spike density function\n",
    "plt.plot(sdf_centers, convolved_spiketimes)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interspike intervals (ISIs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:21:49.680886Z",
     "start_time": "2017-07-26T18:21:49.546422Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's work with the same neuron as before, this time with the spikes from the entire session. \n",
    "# Plot the spikes\n",
    "plt.plot(these_spikes.time, np.ones(len(these_spikes.time)), '|', ms=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:21:57.102338Z",
     "start_time": "2017-07-26T18:21:57.098835Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the duration of the interspike intervals\n",
    "isi = np.diff(these_spikes.time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:21:58.076295Z",
     "start_time": "2017-07-26T18:21:57.827827Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the binned ISIs\n",
    "plt.hist(isi, 50)\n",
    "plt.ylim(0, 120)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spike autocorrelation function (ACF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:22:05.601712Z",
     "start_time": "2017-07-26T18:22:05.579197Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def autocorrelation(spiketimes, bin_size, max_time):\n",
    "    \"\"\"Computes the autocorrelation for an individual spiketrain.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    spiketimes : np.array\n",
    "    bin_size : float\n",
    "    max_time : float\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    aurocorrelation : np.array\n",
    "    bin_centers : np.array    \n",
    "    \"\"\"\n",
    "    \n",
    "    bin_centers = np.arange(-max_time-bin_size, max_time+bin_size, bin_size)\n",
    "    autocorrelation = np.zeros(bin_centers.shape[0]-1)\n",
    "    \n",
    "    for spike in spiketimes:\n",
    "        relative_spike_time = spiketimes - spike\n",
    "        \n",
    "        autocorrelation += np.histogram(relative_spike_time, bin_centers)[0]\n",
    "    \n",
    "    bin_centers = bin_centers[2:-1]\n",
    "    autocorrelation = autocorrelation[1:-1]\n",
    "    \n",
    "    # Normalize the autocorrelation\n",
    "    autocorrelation /= np.max(autocorrelation)\n",
    "        \n",
    "    return autocorrelation, bin_centers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:22:08.024446Z",
     "start_time": "2017-07-26T18:22:07.820764Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the autocorrelation for our spikes of interest\n",
    "acf, bin_centers = autocorrelation(these_spikes.time, bin_size=0.001, max_time=1.)\n",
    "\n",
    "# Plot the autocorrelation\n",
    "plt.plot(bin_centers, acf)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spike cross-correlation function (CCF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:22:11.490283Z",
     "start_time": "2017-07-26T18:22:11.467758Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def crosscorrelation(spiketimes1, spiketimes2, bin_size, max_time):\n",
    "    \"\"\"Computes the autocorrelation for an individual spiketrain.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    spiketimes1 : np.array\n",
    "    spiketimes2 : np.array\n",
    "    bin_size : float\n",
    "    max_time : float\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    crosscorrelation : np.array\n",
    "    bin_centers : np.array    \n",
    "    \"\"\"\n",
    "    \n",
    "    bin_centers = np.arange(-max_time-bin_size, max_time+bin_size, bin_size)\n",
    "    crosscorrelation = np.zeros(bin_centers.shape[0]-1)\n",
    "    \n",
    "    for spike in spiketimes1:\n",
    "        relative_spike_time = spiketimes2 - spike\n",
    "        \n",
    "        crosscorrelation += np.histogram(relative_spike_time, bin_centers)[0]\n",
    "    \n",
    "    bin_centers = bin_centers[2:-1]\n",
    "    crosscorrelation = crosscorrelation[1:-1]\n",
    "    \n",
    "    # Normalize the crosscorrelation by the number of spikes in the first input\n",
    "    crosscorrelation /= len(spiketimes1)\n",
    "        \n",
    "    return crosscorrelation, bin_centers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:22:12.673015Z",
     "start_time": "2017-07-26T18:22:12.363520Z"
    },
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# limit spikes to task times\n",
    "task_start = info.task_times['task'].start\n",
    "task_stop = info.task_times['task'].stop\n",
    "\n",
    "task_spikes = [spiketrain.time_slice(task_start, task_stop) for spiketrain in spikes]\n",
    "\n",
    "# Find the crosscorrelation for our spikes of interest\n",
    "idx1 = 73\n",
    "idx2 = 31\n",
    "\n",
    "ccf, bin_centers = crosscorrelation(task_spikes[idx1].time, task_spikes[idx2].time, \n",
    "                                    bin_size=0.01, max_time=1.)\n",
    "\n",
    "# Plot the crosscorrelation\n",
    "plt.plot(bin_centers, ccf)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating fake spike data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:22:39.580026Z",
     "start_time": "2017-07-26T18:22:39.418912Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate a fake spike train\n",
    "dt = 0.01\n",
    "spiketime = np.arange(-1, 1, dt)\n",
    "\n",
    "probability = 0.5\n",
    "random_values = np.random.random((1, len(spiketime)))\n",
    "spike_idx = np.where(random_values < probability)[1]\n",
    "toy_spikes = nept.SpikeTrain(spiketime[spike_idx])\n",
    "\n",
    "# Plot the fake spikes\n",
    "plt.plot(toy_spikes.time, np.ones(len(toy_spikes.time)), '|', ms=30, color='k')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "299px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}