septum-mec/actions/stimulus-response/data/10-calculate-stimulus-response.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "07:55:15 [I] klustakwik KlustaKwik2 version 0.2.6\n",
      "/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:20: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
      "Please use `tqdm.notebook.*` instead of `tqdm._tqdm_notebook.*`\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import spatial_maps as sp\n",
    "import septum_mec.analysis.data_processing as dp\n",
    "import septum_mec.analysis.registration\n",
    "import expipe\n",
    "import os\n",
    "import pathlib\n",
    "import numpy as np\n",
    "import exdir\n",
    "import pandas as pd\n",
    "import optogenetics as og\n",
    "import quantities as pq\n",
    "import shutil\n",
    "from distutils.dir_util import copy_tree\n",
    "\n",
    "from septum_mec.analysis.stimulus_response import stimulus_response_latency, compute_response\n",
    "\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "from tqdm._tqdm_notebook import tqdm_notebook\n",
    "tqdm_notebook.pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_gaussian_kde = 0.04\n",
    "window_size = 0.03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_loader = dp.Data()\n",
    "actions = data_loader.actions\n",
    "project = data_loader.project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = pathlib.Path('output/stimulus-response')\n",
    "(output / 'data').mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "identify_neurons = actions['identify-neurons']\n",
    "units = pd.read_csv(identify_neurons.data_path('units'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def process(row):\n",
    "    \n",
    "    action_id = row['action']\n",
    "    channel_id = int(row['channel_group'])\n",
    "    unit_id = int(row['unit_name'])    \n",
    "    \n",
    "    spike_times = data_loader.spike_train(action_id, channel_id, unit_id)\n",
    "    \n",
    "    spike_times = np.array(spike_times)\n",
    "    \n",
    "    stim_times = data_loader.stim_times(action_id)\n",
    "    \n",
    "    nan_series = pd.Series({\n",
    "            't_e_peak': np.nan,\n",
    "            'p_e_peak': np.nan,\n",
    "            't_i_peak': np.nan,\n",
    "            'p_i_peak': np.nan\n",
    "        })\n",
    "    \n",
    "    if stim_times is None:\n",
    "        return nan_series\n",
    "    \n",
    "    stim_times = np.array(stim_times)\n",
    "    \n",
    "    times, spikes, kernel, p_e, p_i = stimulus_response_latency(\n",
    "        spike_times, stim_times, window_size, std_gaussian_kde)\n",
    "    \n",
    "    # if no spikes detected after stimulus nan is returned\n",
    "    if all(np.isnan([p_e, p_i])):\n",
    "        return nan_series\n",
    "        \n",
    "    t_e_peak, p_e_peak, t_i_peak, p_i_peak = compute_response(\n",
    "        spike_times, stim_times, times, kernel, p_e, p_i)\n",
    "\n",
    "    return pd.Series({\n",
    "        't_e_peak': t_e_peak,\n",
    "        'p_e_peak': p_e_peak,\n",
    "        't_i_peak': t_i_peak,\n",
    "        'p_i_peak': p_i_peak\n",
    "    })\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "476c31da67274b2396ed3f2ec54a8344",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1298), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mikkel/apps/expipe-project/septum-mec/septum_mec/analysis/stimulus_response.py:33: RuntimeWarning: invalid value encountered in less\n",
      "  if any(times[idxs_i] < te_peak):\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "results = units.merge(\n",
    "    units.progress_apply(process, axis=1), \n",
    "    left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEFCAYAAADpIfy5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHZlJREFUeJzt3XuUXGWZ7/Fvd5JOgwaUMQMSwcjtZ89REFCDB2PiGiBcjgNHPcIgeDuIc8yoLG9kHJB4QUUDDk4QHDBGGVl6hOCoTCRzUDgxojkqqMH2YUABZ1wuCSMQZdpOQp8/3t2k0nTddlftqtr1+6yVlapdtXc9u/rdT7/97vcyMDExgZmZlc9gpwMwM7P2cII3MyspJ3gzs5JygjczKykneDOzknKCNzMrKSf4PiNpraR3dzoO6y+S7pT0tII/81ZJry7yM7vN7E4HYGblFxEv6HQM/cgJPidJS4FPAP8OHAT8J/CGiBitsc8IcDnwJ8As4FMRsabO56wFJoARYD6wAXh7RGyvdjxJg8AngWOAecAAcE5EbJpy7MuAI4BTI+L3TX0B1rdylv0JYH5EbK3xnluBnwEvBJ4BXBsRF2Wv/VfgEuApwOPAyoj4hqSnAFcChwH7ANuAMyMiKo47G7gO2A68PiJ25Dvz3uMmmpk5Crg0Ig4HPgdcW+2NWSG7HlgREUcDS4B3Szqmgc85AjgO+LPs31vqHG8RsD/wkoj4M+DzwIqK4w1IugJYCJzs5G45NFz2m/Rs4Njs+KdL+m+Snp59xtkRcRTwF8CVkg4ETgIejohjIuIw4P8Bf11xvCHgK8BvgbP6KbmDa/Az9eOI2Jg9XgNcIelPIuKhad57GHAwsEbS5LY9gCOB79X5nLWTSVjSF4DTgG9VO15EXCnpAtIvgoOBpaSazaR3An8KvCAi/tjoyZpVaKbsN+MzEbEdeFjSV4BlpBr7M4GvVpT1CeDwiLhe0i8kvQ04hFTWb6843qWkv2IPjoi+m5fFCX5mKmsDA9m/nVXeO4tU03iiLVLSvsAjTX7OYPYZVY8n6RRS082lwD8BPwfOqjjGbcAmYK2kY7ILyqwZzZT9vMetLOujEbFo8gVJ+wMPSvpfwLnAalIzzH8Az6k4xrVZbFeTav59xU00M/MCSYdnj88FNkXEw1XeG8CYpLMAJB0AbAGObuBzTpc0V9Iw8Hrg63WOdzzw9Yi4kvQn62mki2TSD0gXxMPAygbP1axSM2W/GWdJGsyaZV5DKuvfAw6V9DIASS8A/pXUDLmM9BfuZ0nXxCvYvaxvBi4EDpH05hbE11Oc4GfmN8DFkn5KSqJnV3tjRIwDpwLnSPoJ6WbphVNvfFbxGLAR+Gn2/+fqHO8qYEm2/XbgXuA52c3XyXgmgDcBb81uYJk1o+Gy36Q9SEn5e8CnI+KWiHgQeBXwCUk/JtXKz46I+4FVpKbIO4FbgB+RmmqeEBFjwBuy/Q9uUZw9YcDTBeeT9SRYHRHPa/PnrAW2RMSqdn6OWaPaVfazXjSrI+L6Vh63n7kNvoUkvQd4bZWXPxERX5xmHwFfrrJPkLqgmXW1Rsq+pJeTuu9O59vtiay/uQZvZlZSboM3MyuphppoJP0p8ENS74wdwFpSP9QtwPKIeFzSRcAp2evnRcTmtkRsZmYNqZvgJc0BPsOutuDLgAsi4lZJVwGnSrqfNJJyEXAAcAPwonrHXrRo0cSCBQuetH18fJyhoaGGT6LM/F3sMt13cdddd22NiPkdCqlp1co89MbPuhdihN6IcyYxNlruG6nBryJ1u/ub7PnRpIEyAOuBE0g3AzdkXe8ekDRb0vyse1NV8+fP5+KLL37S9rGxMYaHhxsIrfz8Xewy3Xdx2mmn3d+hcHJZsGAB69atm/a10dFRRkZGCo6oOb0QI/RGnDOJMatU11UzwUt6A/BgRNwsaTLBD1QM+d0G7A3sBVQOUZ7cXjPBDw0NTXuCvfDDKYq/i138XZg1p14N/k3AhKTjgBcAXyDNYTJpHmk05KPZ46nbzcysQ2r2oomIl0XEkohYCtwJvA5Ynw10gDST20bSvCbLsiHGBwKDtaYFNTOz9ssz0OldwNWShoBR4PqI2ClpI2lY/CCwfCZBHbjwoFz7jW3fyfCcWfXfaNaF8pR7l3mrpeEEn9XiJy2Z5vWVtGjiqqfsMZeFK25qer/7PnZKKz7erCPylHuXeavFA53MzErKCd7MrKSc4M3MSsoJ3syspJzgzcxKygnezKyknODNzErKKzqZNSibm+kN2dNh0vQdf0makO9X2faLIuK2J+1s1gFO8GYNioi1pLUQkHQFsIY0u+p7I+KGzkVmNj0neLMmSXoh8F8iYrmk9cCRks4DNgPnR8SOavuOj48zOjo67Wt5Z8qsdrx2GBsbK/Tz8uqFOIuI0QnerHnvAz6QPf4X4KvAL0nrJvwVsLrajtWmyJ6JIqdQ7pUpm3shziJi9E1WsyZIehqgiPh2tmlNRPwiWyPhn4AjOxed2e6c4M2a8zLgFgBJA8BPJD0re+3PSWsXm3UFJ3iz5gj4BUBWaz8HWCfpNmBP4OoOxma2G7fBmzUhIj4x5fkGYEOHwjGryTV4M7OScoI3MyspJ3gzs5JygjczKykneDOzknKCNzMrKSd4M7OScoI3MyspJ3gzs5JygjczKykneDOzknKCNzMrKSd4M7OScoI3MyspJ3gzs5JygjfrYWPbdxayj/UmL/hh1sOG58xi4Yqbmtrnvo+d0qZorNu4Bm9mVlKuwZs1SNKPgEezp78EPgNcDuwANkTEBzoVm9l0nODNGiBpGBiIiKUV2+4EXkVahPsmSUdGxB0dCtHsSZzgzRpzBLCnpA2k62YlMDci7gWQdDNwHOAEb13DCd6sMY8Bq4BrgEOB9cDDFa9vAw6qd5Dx8XFGR0enfW1kZGTmUTaoWgz1jI2N5d63SL0QZxExOsGbNeZu4J6ImADulvQIsE/F6/PYPeFPa2hoqNBEXk3eGEZHR7si/np6Ic4iYqyZ4CXNAdYAC4G5wIeBnwFrgQlgC7A8Ih6XdBFwCumG03kRsbl9YZsV7k3A84G3Stof2BP4g6SDSW3wywDfZLWuUq+b5FnAQxGxGDgRWA1cBlyQbRsATpV0FLAEWAScAVzRvpDNOuKzwNMkfQf4MinhnwN8EdgM3BER3+9gfGZPUq+J5ivA9dnjAVLt/GjgtmzbeuAEIEjdxCaAByTNljQ/Ih5sQ8xmhYuIceDMaV46puhYzBpVM8FHxO8BJM0jJfoLgFVZIod0Y2lvYC/goYpdJ7fXTPDVbjjNpF2q22+sNKsXbhYVxd+FWXPq3mSVdABwI/DpiLhO0scrXp68sfRo9njq9praccOp22+sNKsXbhYVxd+FWXNqtsFL2hfYAJwfEWuyzXdIWpo9PgnYCGwClkkalHQgMBgRW9sUs5mZNaBeDf59wNOBCyVdmG17B/ApSUPAKHB9ROyUtBG4nfRLY3m7AjYzs8bUa4N/BymhT7VkmveuJI3u65ix7TsZnjOr7fuYmfWCUg108tSpZma7eLpgM7OScoI3MyspJ3gzs5JygjczKykneDOzknKCNzMrKSd4M7OScoI3MyspJ3gzs5JygjczKykneDOzknKCNzMrqVJNNmbWLlUWoP8V8A3gX7O3XRkRX+5IgGbTcII3a8zkAvRnS9oHuBP4IHBZRFza2dDMpucEb9aYagvQS9KppFr8eRGxrUPxNczrJvQPJ3izBlRZgH4ucE1E/FDS3wIXAe+udZxqC81DcesJ5103YXR0tGcWPu+FOIuI0QnerEHTLED/tIiYXFz+RuDv6x2jHQvNF2VkZKRnFj7vhTiLiNG9aMwaUGUB+pslvTh7/OfADzsSnFkVrsGbNWa6BejfCXxS0nbgN8C5nQrObDpO8GYNqLEA/bFFx2LWKDfRmJmVlBO8mVlJOcGbmZWUE7yZWUk5wZtZXWPbdwLND8aa3M86w71ozKyuPKNfIY2Atc7p+xp8nhqGayVm1gv6vgafd14OM6vPE5t1Vt8neDNrH1egOqvvm2jMzMrKCd7
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results.loc[:, ['t_e_peak', 't_i_peak', 'p_e_peak', 'p_i_peak']].hist()\n",
    "plt.gcf().savefig(output / 'figures' / 'summary_histogram.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:2: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
      "Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2828a0134cc2449ba15a30119991324d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1284), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "psth, time = {}, {}\n",
    "for i, row in tqdm(units.iterrows(), total=len(units)):\n",
    "    action_id = row['action']\n",
    "    channel_group = row['channel_group']\n",
    "    unit_name = row['unit_name']\n",
    "    name = f'{action_id}_{channel_group}_{unit_name}'\n",
    "    spike_times = data_loader.spike_train(action_id, channel_group, unit_name)\n",
    "    \n",
    "    spike_times = np.array(spike_times)\n",
    "    \n",
    "    stim_times = data_loader.stim_times(action_id)\n",
    "    \n",
    "    if stim_times is None:\n",
    "        continue\n",
    "    \n",
    "    stim_times = np.array(stim_times)\n",
    "    \n",
    "    _, spikes, kernel, p_e, p_i = stimulus_response_latency(\n",
    "        spike_times, stim_times, window_size, std_gaussian_kde)\n",
    "    \n",
    "    if all(np.isnan([p_e, p_i])):\n",
    "        continue\n",
    "    times = np.arange(-0.005, window_size, 1e-4)\n",
    "    psth.update({name: kernel(times)})\n",
    "    time.update({name: times})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(psth).to_feather(output / 'data' / 'psth.feather')\n",
    "pd.DataFrame(time).to_feather(output / 'data' / 'times.feather')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save to expipe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "action = project.require_action(\"stimulus-response\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "action.modules['parameters'] = {\n",
    "    'window_size': window_size,\n",
    "    'std_gaussian_kde': std_gaussian_kde\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "action.data['results'] = 'results.csv'\n",
    "results.to_csv(action.data_path('results'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/media/storage/expipe/septum-mec/actions/stimulus-response/data/data/times.feather',\n",
       " '/media/storage/expipe/septum-mec/actions/stimulus-response/data/data/psth.feather',\n",
       " '/media/storage/expipe/septum-mec/actions/stimulus-response/data/figures/summary_histogram.png']"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "copy_tree(output, str(action.data_path()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "septum_mec.analysis.registration.store_notebook(action, \"10-calculate-stimulus-response.ipynb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}