{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 05 - Tuning curves and decoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals\n",
    "\n",
    "- Learn to estimate and plot 2D tuning curves\n",
    "- Implement a Bayesian decoding algorithm\n",
    "- Compare the decoded and actual positions by computing the decoding error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the tuning curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:19.560027Z",
     "start_time": "2017-07-27T14:50:18.594812Z"
    },
    "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",
    "\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-27T14:50:19.571017Z",
     "start_time": "2017-07-27T14:50:19.562010Z"
    },
    "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-27T14:50:20.110400Z",
     "start_time": "2017-07-27T14:50:19.573019Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load position (.nvt) from this experiment\n",
    "position = nept.load_position(os.path.join(data_folder, info.position_filename), info.pxl_to_cm)\n",
    "\n",
    "# Plot the position\n",
    "plt.plot(position.x, position.y, 'k.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:20.784880Z",
     "start_time": "2017-07-27T14:50:20.111902Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load spikes (.t and ._t) from this experiment\n",
    "spikes = nept.load_spikes(data_folder)\n",
    "\n",
    "# Plot the spikes\n",
    "for idx, spiketrain in enumerate(spikes):\n",
    "    plt.plot(spiketrain.time, np.ones(len(spiketrain.time))+idx, '|')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:20.802893Z",
     "start_time": "2017-07-27T14:50:20.786882Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# limit position and spikes to task times\n",
    "task_start = info.task_times['task'].start\n",
    "task_stop = info.task_times['task'].stop\n",
    "\n",
    "task_position = position.time_slice(task_start, task_stop)\n",
    "task_spikes = [spiketrain.time_slice(task_start, task_stop) for spiketrain in spikes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:20.815902Z",
     "start_time": "2017-07-27T14:50:20.804895Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# limit position to those where the rat is running\n",
    "run_position = nept.speed_threshold(task_position, speed_limit=1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:20.971013Z",
     "start_time": "2017-07-27T14:50:20.817403Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the running Y position over time\n",
    "plt.plot(run_position.time, run_position.y, 'b.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:21.128143Z",
     "start_time": "2017-07-27T14:50:20.974014Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the running position\n",
    "plt.plot(run_position.x, run_position.y, 'b.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:21.609967Z",
     "start_time": "2017-07-27T14:50:21.129625Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the task spikes\n",
    "for idx, spiketrain in enumerate(task_spikes):\n",
    "    plt.plot(spiketrain.time, np.ones(len(spiketrain.time))+idx, '|', color='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:21.619715Z",
     "start_time": "2017-07-27T14:50:21.611468Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the X and Y boundaries from the unfiltered position, with 3 cm bins\n",
    "xedges, yedges = nept.get_xyedges(position, binsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:23.798547Z",
     "start_time": "2017-07-27T14:50:21.621475Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tuning_curves = nept.tuning_curve_2d(run_position, np.array(task_spikes), xedges, yedges, \n",
    "                                     occupied_thresh=0.2, gaussian_sigma=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:24.497020Z",
     "start_time": "2017-07-27T14:50:23.800025Z"
    },
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot a few of the neuron's tuning curves\n",
    "xx, yy = np.meshgrid(xedges, yedges)\n",
    "\n",
    "for i in [7, 33, 41]:\n",
    "    print('neuron:', i)\n",
    "    plt.figure(figsize=(6, 5))\n",
    "    pp = plt.pcolormesh(xx, yy, tuning_curves[i], cmap='bone_r')\n",
    "    plt.colorbar(pp)\n",
    "    plt.axis('off')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decoding\n",
    "\n",
    "Next, let's decode the location of the subject using a Bayesian algorithm.\n",
    "\n",
    "Specifically, this is a method known as \"one-step Bayesian decoding\" and is\n",
    "illustrated in this figure from van der Meer et al., 2010.\n",
    "\n",
    "<img src=\"05-bayesian-decoding.png\" width=\"85%\" height=\"85%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:25.589313Z",
     "start_time": "2017-07-27T14:50:24.499021Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Bin the spikes\n",
    "window_size = 0.0125\n",
    "window_advance = 0.0125\n",
    "\n",
    "time_edges = nept.get_edges(run_position, window_advance, lastbin=True)\n",
    "counts = nept.bin_spikes(task_spikes, run_position, window_size, window_advance,\n",
    "                         gaussian_std=None, normalized=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:25.601050Z",
     "start_time": "2017-07-27T14:50:25.591298Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reshape the 2D tuning curves (essentially flatten them, while keeping the 2D information intact)\n",
    "tc_shape = tuning_curves.shape\n",
    "decode_tuning_curves = tuning_curves.reshape(tc_shape[0], tc_shape[1] * tc_shape[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:53.900455Z",
     "start_time": "2017-07-27T14:50:25.603306Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the likelihoods\n",
    "likelihood = nept.bayesian_prob(counts, decode_tuning_curves, window_size, min_neurons=2, min_spikes=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:53.909962Z",
     "start_time": "2017-07-27T14:50:53.902957Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the center of the position bins\n",
    "xcenters = (xedges[1:] + xedges[:-1]) / 2.\n",
    "ycenters = (yedges[1:] + yedges[:-1]) / 2.\n",
    "xy_centers = nept.cartesian(xcenters, ycenters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:56.555359Z",
     "start_time": "2017-07-27T14:50:53.912964Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Based on the likelihoods, find the decoded location\n",
    "decoded = nept.decode_location(likelihood, xy_centers, time_edges)\n",
    "nan_idx = np.logical_and(np.isnan(decoded.x), np.isnan(decoded.y))\n",
    "decoded = decoded[~nan_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:56.719959Z",
     "start_time": "2017-07-27T14:50:56.556843Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the decoded position\n",
    "plt.plot(decoded.x, decoded.y, 'r.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the decoded to actual positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:56.727983Z",
     "start_time": "2017-07-27T14:50:56.721460Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the actual position for every decoded time point\n",
    "actual_x = np.interp(decoded.time, run_position.time, run_position.x)\n",
    "actual_y = np.interp(decoded.time, run_position.time, run_position.y)\n",
    "actual_position = nept.Position(np.hstack((actual_x[..., np.newaxis],\n",
    "                                           actual_y[..., np.newaxis])), decoded.time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:56.916600Z",
     "start_time": "2017-07-27T14:50:56.729466Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the actual position\n",
    "plt.plot(actual_position.x, actual_position.y, 'g.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the pedestal is not represented as round, as before. \n",
    "This is because we are interpolating to find an actual position \n",
    "that corresponds to each decoded time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-27T14:50:57.096253Z",
     "start_time": "2017-07-27T14:50:56.918601Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Find the error between actual and decoded positions\n",
    "errors = actual_position.distance(decoded)\n",
    "print('Mean error:', np.mean(errors), 'cm')\n",
    "\n",
    "# Plot the errors\n",
    "plt.hist(errors)\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": "281px",
    "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
}