Science
/
amgen
1
1
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Topology in neuroscience
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
20 Commits
1 Branch
0 Tags
2.3 GiB
Tree: b01ba9aba1
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'b01ba9aba1'
${ noResults }
amgen/exp/10_PSTH-STC.ipynb

774 lines
176 KiB
Raw Normal View History Unescape

Processing experimental data for flashed gratings
2 years ago
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"stim_data.pkl\n",
"dict_keys(['stim_val', 'trial_stim_id', 'key_list', 'num_trial', 'trial_pair_id', 'pair_val', 'pair_trial_id', 'stim_id_trial', 'num_stim'])\n",
"spike_data.pkl\n",
"dict_keys(['spike_count_rate', 'avg_firing_rate', 'sem_firing_rate', 'firing_rate', 'stim_num_trial', 'C_r_fphi_theta', 'theta_hist', 'phase_hist', 'pair_hist'])\n",
"corr_data.pkl\n",
"dict_keys(['corr_stim_unit', 'optimal_avg_firing_rate', 'stim_hist', 'stim_hist_caution'])\n"
]
}
],
"source": [
"# load data\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import scipy.stats as sts\n",
"import pickle\n",
"import ipywidgets as widgets\n",
"from mpl_toolkits import mplot3d\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from scipy import sparse\n",
"from matplotlib import cm\n",
"cmap_hot = cm.get_cmap('hot')\n",
"cmap_viridis = cm.get_cmap('viridis')\n",
"\n",
"execfile('Stimulus.py')\n",
"data_folder = \"../data/\"\n",
"\n",
"stim_file = \"Stiminfo_PVCre_2021_0012_s06_e14.csv\"\n",
"stim = pd.read_csv(data_folder+stim_file)\n",
"\n",
"spike_times_file = \"Spiketimes_PVCre_2021_0012_s06_e14.npy\"\n",
"spike_times = np.load(data_folder+spike_times_file, allow_pickle=True)\n",
"active = [len(spike_times[i]) > 0 for i in range(len(spike_times))]\n",
"spike_times = spike_times[np.where(active)]\n",
"\n",
"num_unit = len(spike_times)\n",
"num_trial = len(stim)\n",
"\n",
"# sort by firing rate\n",
"\n",
"num_spike = list(map(len, spike_times))\n",
"# num_spike = np.array([len(spike_times[i]) for i in range(len(spike_times))])\n",
"spike_times = spike_times[np.argsort(num_spike)[::-1]]\n",
"execfile('load.py')\n",
"\n",
"max_delay = 300 # dt\n",
"tau_id_range = np.arange(max_delay)\n",
"\n",
"latest_spike_time = max([np.max(s) for s in spike_times if len(s)])\n",
"latest_stim_offtime = list(stim['stim_offtime'])[-1]\n",
"experiment_dur = max([latest_spike_time, latest_stim_offtime])\n",
"\n",
"dt = 0.001 # 1 ms\n",
"exp_time = np.arange(0, experiment_dur, dt)\n",
"trial_length = np.mean(stim['stim_offtime']-stim['stim_ontime'])\n",
"M = len(exp_time)\n",
"\n",
"# binary spike and stimulus trains\n",
"B_stim = {}\n",
"for key in key_list:\n",
" B_stim[key] = []\n",
" for stim_id, trials in enumerate(stim_id_trial[key]):\n",
" B_stim[key].append([])\n",
" s = []\n",
" for trial_id in trials:\n",
" t_on, t_off = stim['stim_ontime'][trial_id], stim['stim_offtime'][trial_id]\n",
" s += list(np.arange(int(t_on//dt), int(t_off//dt)))\n",
"\n",
" B_stim[key][stim_id] = sparse.coo_matrix((np.ones(len(s)), (np.zeros(len(s), dtype=int), s)), shape=(1, M))\n",
"s = spike_times//dt\n",
"B_spike = []\n",
"for unit_id in range(num_unit):\n",
" B_spike.append(sparse.coo_matrix((np.ones(len(s[unit_id])), (np.zeros(len(s[unit_id]), dtype=int), np.int0(s[unit_id]))), shape=(1, M)))\n",
"\n",
"# histogram error bars: num spikes\n",
"s = np.zeros((num_unit, 2))\n",
"for unit_id in range(num_unit):\n",
" # print(\"unit: %d\"%unit_id)\n",
" a = np.zeros(len(tau_id_range))\n",
" for tau_id in tau_id_range:\n",
" a[tau_id] = np.sum(B_spike[unit_id].col >= tau_id)\n",
" \n",
" s[unit_id] = [np.mean(a), np.std(a)]\n",
"\n",
"key_symbol = {'pair':'$(\\\\theta,\\phi)$', 'orientation':'$\\\\theta$', 'phase':'$\\phi$'}\n",
"\n",
"# 2D tuning\n",
"avg_firing_rate_pair = np.array([sts.zscore(stim_hist['pair'][unit_id]).reshape((len(tau_id_range), num_stim['orientation'], num_stim['phase'])) for unit_id in range(num_unit)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Peristimulus time histogram\n",
"we should wait for the PSTH to stabilize otherwise we get the transient sharp response to sudden change of stimulus."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# concatenate spike times of all trials with the same stimuli\n",
"# firing rate"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/lib/python3.10/site-packages/pandas/core/roperator.py:13: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
" return right - left\n"
]
}
],
"source": [
"trial_spike_times = np.ndarray((num_unit, num_trial), dtype=object)\n",
"delta_t_on = 0.2 # 200ms\n",
"delta_t_off = 0.3\n",
"for trial_id in range(num_trial):\n",
" t_on, t_off = stim['stim_ontime'][trial_id]-delta_t_on, stim['stim_offtime'][trial_id]+delta_t_off\n",
" for unit_id in range(num_unit):\n",
" trial_spike_times[unit_id, trial_id] = spike_times[unit_id][np.where((spike_times[unit_id] < t_off) & (spike_times[unit_id] > t_on))]\n",
"subseq_trials = 5\n",
"subseq_spike_times = np.ndarray(num_unit, dtype=object)\n",
"for unit_id in range(num_unit):\n",
" subseq_spike_times[unit_id] = [np.concatenate(trial_spike_times[unit_id][i:i+subseq_trials]) for i in range(num_trial-subseq_trials)]\n",
"\n",
"subseq_spike_times_locked = {}\n",
"subseq_spike_times_locked['ontime'] = [subseq_spike_times[unit_id] - stim['stim_ontime'][:-subseq_trials] for unit_id in range(num_unit)]\n",
"subseq_spike_times_locked['offtime'] = [subseq_spike_times[unit_id] - stim['stim_offtime'][:-subseq_trials] for unit_id in range(num_unit)]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import zscore\n",
"trial_dur = stim['stim_offtime'] - stim['stim_ontime']\n",
"z_trial_dur = zscore(trial_dur)\n",
"intertrial_dur = np.array(stim['stim_ontime'])[1:] - np.array(stim['stim_offtime'])[:-1]\n",
"z_intertrial_dur = zscore(intertrial_dur)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAANSCAYAAABIrSnGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABvY0lEQVR4nO3de5wkaV3n++8vIjOrqqt7bl09N+ZWI8hVHI4to4J7xgsIHEUQXr7APYjHC3h0jrqr5zi6e2RcebHIgrhn2WUXhOVyFLxyRIVVvKLuvththtsMoMJ0j8wwzHR1z62ruiszMp7zR0RkRkRGZGZFVWVUZnzer1e/KjPyiSeeqoKa7/PEE89jzjkBAAAA2Bmv7gYAAAAA84ggDQAAAFRAkAYAAAAqIEgDAAAAFRCkAQAAgAoI0gAAAEAFUwVpM3unmT1oZnemjt1uZveZ2Sfjfy/InXOdmZ0zs59JHTtlZp+Jy59IHb/MzD5iZv8Qf710L745AAAAYL9MOyL9LknPKzj+ZufcTfG/D+U++xVJHy4451vi8sdTx26T9GfOuSdI+rP4PQAAAHBgTRWknXMflXR22krN7EWSTkq6a8pTvlvSu+PX75b0ommvBQAAANShtcvzbzWz75d0QtJPO+ceMrPDkn5W0nMk/UyuvJP0J2bmJP0n59zb4uNXOOfuj19/RdIVRRczs1dJepUkra6uft2TnvSkXTYfAJrlzvse0drhJV158XLdTQGAufHxj398wzl3LH98N0H6rZJ+SVE4/iVJb5L0g5JuVzTl45yZ5c95tnPuPjO7XNJHzOzz8Wj3gHPOxUF7RBy83yZJx48fdydOnCgqBgAo8cR/+WH9wLNu0M89/8l1NwUA5oaZ3VN0vHKQds49kKr87ZL+MH57s6SXmtkbJF0iKTSzC865tzjn7ovPfdDMPiDpmZI+KukBM7vKOXe/mV0l6cGq7QIAAABmofLyd3HgTbxY0p2S5Jz7ZufcDc65GyT9qqTXOefeYmarZnYkPndV0nOTcyR9UNIr49evlPT7VdsFAAAAzMJUI9Jm9j5Jt0haM7N7Jb1G0i1mdpOiqR2nJL16QjVXSPpAPN2jJek3nHP/Jf7s9ZJ+y8x+SNI9kr53R98FAAAAMGNTBWnn3MsLDr9jivNuT72+W9LXlpQ7I+nbpmkLAAAAcBCwsyEAAABQAUEaAAAAqIAgDQAAAFRAkAYAAAAqIEgDAAAAFRCkAQAAgAoI0gAAAEAFBGkAAACgAoI0AAAAUAFBGgAAAKiAIA0ADWImhaGruxkAsBAI0gDQIFdctKwvP3Kh7mYAwEIgSANAg9y4tqqTpzfrbgYALASCNAA0yPraYZ3c2JRzTO8AgN0iSANAg6yvHdL5Xl8PPrZdd1MAYO4RpAGgQQ4vtyRJ57v9mlsCAPOPIA0AAABUQJAGAAAAKiBIAwAAABUQpAEAAIAKCNIAAABABQRpAAAAoAKCNAAAAFABQRoAAACogCANAAAAVECQBgAAACogSAMAAAAVEKQBAACACgjSAAAAQAUEaQAAAKACgjQAAABQAUEaABrk4a2eJOnwcqvmlgDA/CNIA0CDnNzY1JHllo6udupuCgDMPYI0ADTIyY1N3bi2KjOruykAMPcI0gDQIHef3tT62mrdzQCAhUCQBoAG2eoGWukwPxoA9gJBGgAa5Ia1Vd1zZrPuZgDAQiBIA0CDrK+t6uQGQRoA9gJBGgAa5Ma1Vd3/yAVtdYO6mwIAc48gDQANctXFK5KkBx/drrklADD/CNIA0CAef/UBYM/wJxUAAACogCANAAAAVECQBoAG2er2JUmdFn/+AWC3+EsKAA1yz5ktdVqerrxoue6mAMDcI0gDQIPcfXpT60dX5XlWd1MAYO4RpAGgQU5unNP62mrdzQCAhUCQBoAG+fLDF3TNpSt1NwMAFgJBGgAaJHROvs+0DgDYCwRpAAAAoAKCNAAAAFABQRoAAACogCANAAAAVECQBgAAACogSAMAAAAVEKQBAACACgjSAAAAQAUEaQAAAKACgjQAAABQAUEaAAAAqIAgDQAAAFRAkAYAAAAqIEgDAAAAFRCkAQAAgAoI0gAAAEAFBGkAAACgAoI0AAAAUAFBGgAAAKiAIA0AAABUQJAGAAAAKiBIAwAAABUQpAEAAIAKCNIAAABABQRpAAAAoAKCNAAAAFABQRoAAACogCANAAAAVECQBgAAACogSAMAAAAVEKQBAACACgjSAAAAQAUEaQAAAKACgjQAAABQAUEaAAAAqIAgDQAAAFRAkAYAAAAqIEgDAAAAFRCkAQAAgAoI0gDQIIc6vh67ENTdDABYCARpAGiQG9ZWdWpjs+5mAMBCIEgDQIOsr63qJEEaAPYEQRoAGuTGtVXd/8gFbXWZ3gEAu0WQBoAGuWilLUnMkwaAPTAxSJvZO83sQTO7M3XsdjO7z8w+Gf97Qe6c68zsnJn9TPz+WjP7CzP7rJndZWY/OW1dAIC9c8+ZLa20fR07vFR3UwBg7rWmKPMuSW+R9J7c8Tc7595Ycs6vSPpw6n0g6aedc3eY2RFJHzezjzjnPjtFXQCAPXJyY1M3rK3K86zupgDA3Js4Iu2c+6iks9NWaGYvknRS0l2pOu53zt0Rv35M0uckPW6njQUA7M7JjU3duLZadzMAYCHsZo70rWb26Xjqx6WSZGaHJf2spF8sO8nMbpD0DEkfG1dXybmvMrMTZnbi9OnTu2g6ADTTmXPbunS1XXczAGAhVA3Sb5X0VZJuknS/pDfFx29XNE3jXNFJcdD+XUk/5Zx7dEJdI5xzb3POHXfOHT927FjFpgNAc1139JC+dPZ83c0AgIUwzRzpEc65B5LXZvZ2SX8Yv71Z0kvN7A2SLpEUmtkF59xbzKytKET/unPu96aoCwCwx9bXDutTX3q47mYAwEKoNCJtZlel3r5Y0p2S5Jz7ZufcDc65GyT9qqTXxSHaJL1D0uecc78yTV0AgL23vraqex/a0nbQr7spADD3Jo5Im9n7JN0iac3M7pX0Gkm3mNlNkpykU5JePaGaZ0l6haTPmNkn42M/75z7kKQ37LAuAEBF11y6otBJDzyyreuOHqq7OQAw1yYGaefcywsOv2OK825Pvf4bSYVrLTnnXjGpLgDA3mj70Z/i0LmaWwIA84+dDQEAAIAKCNIAAABABQRpAGiQbhBKknx2NgSAXSNIA0CD3HNmSy3PdOXFy3U3BQDmHkEaABrk5MamrrvskNo+f/4BYLf4SwoADXJyY1Pra6t1NwMAFgJBGgAa5L6Hz+vqS1bqbgYALASCNAA0yBUXLevBxy7U3QwAWAgEaQBokPW1Vd19erPuZgDAQiBIA0CD3Li2qnvObKkfsrMhAOwWQRoAGuSKi5bV7YfaOLddd1MAYO4RpAGgQe74x4e0drijY4eX6m4KAMw9gjQANESvH+qv/v60vuWJl8tjZ0MA2DWCNAA0xIlTD+mxC4G+7cmX190UAFgIBGkAaIj/+sUNeSY9+wnH6m4KACwEgjQANMShTkuhk5xjxQ4A2AsEaQBoiGRr8JMbrCMNAHuBIA0ADXHjMYI0AOwlgjQANMR1lx2SGUEaAPYKQRoAGmK57avjezrf69fdFABYCARpAAAAoAKCNAAAAFABQRoAAACogCANAAAAVECQBgAAACogSAMAAAAVEKQBAACACgjSAAAAQAUEaQAAAKACgjQAAABQAUEaAAAAqIAgDQAAAFRAkAYAAAAqIEgDAAAAFRCkAQAAgAoI0gAAAEAFBGkAAACgAoI0AAAAUAFBGgAAAKiAIA0AAABUQJAGAAAAKiBIAwAAABUQpAEAAIAKCNIAAABABQRpAAAAoAKCNAAAAFABQRoAAACogCANAAAAVECQBgAAACogSAMAAAAVEKQBAACACgjSAAAAQAUEaQAAAKACgjQAAABQAUEaAAAAqIAgDQAAAFRAkAYAAAAqIEgDAAAAFbTqbgAAYG/0Q6deP1QQOgX9UN1+qKDvFPSdemH02rm6WwkAi4MgDaDxnHPq9Z2CMFSvH4fRfjaUDo7HZdLhNBNaB59H5w7DbKhevq5
"text/plain": [
"<Figure size 864x1080 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(12,15))\n",
"plt.eventplot(intertrial_dur, lineoffsets=-1000, linelengths=1000, color='r')\n",
"plt.plot(intertrial_dur, np.arange(19999))\n",
"plt.ylim([15250, 15450])\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4kAAAJcCAYAAABUquF9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABh6ElEQVR4nO3deXwV5fXH8e9hR8AVBBUxqICKC2pccEHcUax2Fe3iUqvS1v4qXWNt61Y1XWx/1doW21qVn2ttrdZoELe6AQqKSNQokiggASMIGEC28/tjJtebkDA3kJlJmM/bV16589xZzszFm3vueZ5nzN0FAAAAAIAkdUg7AAAAAABA20GSCAAAAADIIUkEAAAAAOSQJAIAAAAAckgSAQAAAAA5JIkAAAAAgJxOaQeQtN69e3tRUVHaYQAAMqq6ulr8HQIApGn69Om17t6nueczlyQWFRVp2rRpaYcBAMio4uJi/g4BAFJlZu9u7Hm6mwIAAAAAckgSAQAAAAA5JIkAAAAAgBySRAAAAABADkkiAAAAACCHJBEAAAAAkEOSCAAAAADIIUkEAAAAAOSQJAIAAAAAckgSAQAAAAA5JIkAAAAAgBySRAAAAABADkkiAAAAACCHJBEAAAAAkEOSCAAAAADIIUkEAAAAAOSQJAIAAAAAckgSAQAAAAA5JIkAAAAAgJx2nySa2SgzqzSz2WZWknY8AAAAANCetesk0cw6SrpZ0imS9pF0tpntk25UAAAAANB+dUo7gM10qKTZ7j5HkszsHklnSHo91aiwxVi33vXJ2nWSJJPl2s0arle/XL/Op8v1z1uD5QbrNN4ZAAAAkKL2niTuImlu3vI8SYelFAu2QLPmL9UZNz+f6DFblGCq4crWzPP5eWiz+23uuI3aT9lvJ133uf0KPyEAAAC0K+09SSyImV0k6SJJGjBgQMrRoD3ZaZtuuuyUveR5bR4ueNjq3nAbd2+0XtPbNbVOfUNz2xSyT22wTmHbNj6PBtuGy09VLtL06iUbrggAAIDYjBk/WZJ078XDEzlee08S50vaNW+5f9jWgLvfIukWSSouLm7iozDQtB237qaLj9kj7TDajLETpquqti7tMAAg05L+sAggfUn//97ek8SXJA0ys4EKksOzJH053ZCALY+76+4X56q8okaSVFRSpnsuOlyH775DypEBQPaQHGYXXxAgKe06SXT3tWZ2iaSJkjpKutXdK1IOC9jiXHjHdD3+xsIGbWfdMkXVpaNTiggAsotEIbt4zZGUdp0kSpK7PyLpkbTjALZkXzlswAZJ4nM/PjalaAAg20gUsosvCJCUdp8kAojfsXvtqOrS0bkxiRPHjUg7JAAAMofkEEnpkHYAAAAAAKKNGT85V00E4kQlEUCkpSvW6Mt/naKK95dJkv72XJUuOGpgylEBAJAtVBKRFCqJACL965V5uQRRkq55+PUUowGAbKOaBCBuVBIBRDrviCLtvdPW+sE/XtW8JSs1+9pT0g4JADKLahKAuJEkAohkZjp89x20787bqEeXTurUkU4IAJAWZrgEEDeSRAAFcXfN/uBjzV70sZ6qXKSj9uytziSLAJA4kkMAcSNJBFCQgZd9ejvS8//+kiSpunR0WuEAAJA5VJGRFMoAAAryzZF7NFj+1Rf2TykSAAAAxIlKIoCC/HjUXqr6oE5VtXWaOG5E2uEAAAAgJiSJAAAA7QhdDrOL1xxJobspgEhzF69QUUmZyitqVLlwuYpKyvTgjPlphwUAAIAYkCQCiHTOrS9u0OaeQiAAAACIHd1NAUR6/HvH6IbHKvXHp9+RJM288iRt3a1zylEBQDbR5TC76GqMpJAkAojUsYPpR6P20pxw4hoSRAAAgC0XSSKAgqxdt17lFTWSpCV1q7Vdjy4pRwQAAIA4MCYRQEGOKH0y9/jAayalGAkAAADiRCURQEHGf+1gfe6PL0iS/uf4QSlHAwBA9jAWEUkhSQRQkAMHbKdRQ/upqrZO3ztxcNrhAEBmMXkJgLiRJAIAALQjJIcA4saYRACR5n+0UkUlZSqvqFHlwuW6eMK0tEMCACBzxoyfnKskA3Gikggg0nfuernB8sSKhSlFAgBAdlFFRlKoJAKI9KevHqxdtu2eW777wsNTjAYAgGyikoikUEkEEKnv1t30fMlxGjthuqpq6zR8jx3SDgkAgMyhkoikUEkEUJAPln+SG5O435UTtaRuddohAUAmUU0CEDcqiQAKcsi1j+ceL1+1VgdeM0nVpaNTjAgAsolqUnZx+xMkhUoigII8+O0jGyw//r0RKUUCAACAOFFJBFCQA3bdVqOG9lNVbZ0mjiNBBAAA2FJRSQQAAAAA5FBJBFCQvz1XpfKKGklSUUmZfn/WMJ0xbJeUowKA7GFcWnbxmiMpVBIBFGTKnA/TDgEAAAAJoJIIoCB/OadYYydMV+XC5XrqByPTDgcAAAAxIUkE0CJdOtIBAQDSRJfD7KKrMZLCpz0ABbn7xfdUXlGjyoXLVVRSpodefT/tkAAAABADkkQAkc659UVd9q/XGrT9z92vpBQNAAAA4kSSCCDSN44a2GD5jGE7a851p6YUDQAAAOLEmEQAkUYM7qPq0tEaO2G6qmrr9PuzDkw7JAAAMoexiEgKlUQAAIB2ZMz4ybkJTAAgDlQSAUR6Y8EynfL7Z3PLRSVlOmHvvvrrucUpRgUA2UQ1CUDcqCQCiPTuhys2aFu7fn0KkQAAkF1UkZEUKokAIo3at1+DMYkTx41IOyQAADKHKjKSQiURAACgHaGaBCBuVBIBRFqzbr1+M7FS5RU1koIxif/85hE6eLftUo4MALKHahKAuFFJBBDp789Xafwzcxq0feFPL6QUDQBkG5VEAHGjkggg0rlHFOn9j1bptheqJUn77LS17hvLN9kAkAYqidlV/+UA/wYQN5JEAJG6duqoK08fqpqlq1RVW6dHvnt02iEBAAAgJnQ3BVAQd9ebNctUuXC5Zs1fqvXrPe2QAAAAEAOSRACRnnu7VgMve0TV4f0ST7vpOV398OspRwUAAIA40N0UQKRunTf8PumCowamEAkAgHFp2cVrjqSQJAKIVFy0vapLR2vshOmqqq3TxHEj0g4JADKLRCG7+IIASSFJRLtz+h+e08x5S9MOI9OKSsrSDgFIXHXp6LRDAAAgEYxJRLtDgggAAADEh0riJnrvwxW6/+V5m7y9RTy/sXkjo7Ztdru8Db2FE1NaCw+av/+Wbpvbrpkz3aFHF31Yt3rTdgoAm2DX7bunHQIAAIlJPEk0s10l3SGpr4Jc6BZ3/72ZXSnpQkkfhKv+xN0fCbe5TNIFktZJ+h93nxi2j5L0e0kdJf3V3UuTOo/P/OE5LV25JqnDAQBSNHfxyrRDAADGIiIxaVQS10r6vru/bGa9JE03s0nhc79z99/kr2xm+0g6S9JQSTtLetzMBodP3yzpREnzJL1kZg+5eyLz8v/nkqP0r1c2rZJYaBWvqQpcSyuAue2a2v9mbFsI24xtN3aiFe8v0xNvLtrUPQPt1p479tTA3j2afG5T3xvCrcPfLS37t/Sg+fvf9ICbO9eoXgubc42+VNx/0zcGWhmTlwCIW+JJorsvkLQgfLzczN6QtMtGNjlD0j3u/omkKjObLenQ8LnZ7j5HkszsnnDdRJLEATtspUtPGBy9IlrVhCnv6oknZ6cdBpCK2Ys+1uxFH+vCowfq8tH7pB0OgJSQHAKIW6oT15hZkaQDJU0Nmy4xs5lmdquZbRe27SJpbt5m88K25tqbOs5FZjbNzKZ98MEHTa2CdqJvr65phwCkboee/H8AZNmY8ZNz1UQAiENqE9eYWU9J/5R0qbsvM7M/SbpGQR+kayTdIOnrrXEsd79F0i2SVFxcvFmdspCuk4b2Yxr6FHGfRABIH5VEAHFLJUk0s84KEsQ73f1fkuTuC/Oe/4ukh8PF+ZJ2zdu8f9imjbQDAAAAWxTGoyIpacxuapL+JukNd/9tXvtO4XhFSfqcpFnh44ck3WVmv1Uwcc0gSS8qmAFhkJkNVJAcniXpy8mcBZAtH61YrSNLn1Td6nWSpKKSMt1z0eE6fPc
"text/plain": [
"<Figure size 1080x720 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"unit_id = 1\n",
"\n",
"fig, ax = plt.subplots(1,2, figsize=(15,10))\n",
"ax[0].plot(np.diff(stim['stim_ontime'], n=2), np.arange(20000-2))\n",
"ax[0].invert_yaxis()\n",
"ax[0].set_xlabel('$\\Delta^2[t_{on}]=\\Delta$ trial duration / s', size=20)\n",
"ax[0].set(ylabel='trials')\n",
"ax[1].eventplot([[trial_dur[i]] for i in range(num_trial)], lineoffsets=1, linelengths=0.8)\n",
"ax[1].set(yticks=[], xlabel='trial duration / s')\n",
"ax[1].invert_yaxis()\n",
"# ax[2].eventplot([[stim['stim_ontime'][i] - stim['stim_offtime'][i]] for i in range(num_trial)], lineoffsets=1, linelengths=0.8)\n",
"# ax[2].set(title='unit %d'%unit_id, yticks=[], xlabel='t/s (locked to trial onsets)')\n",
"# ax[2].invert_yaxis()\n",
"plt.subplots_adjust(wspace=0)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABJ8AAAJsCAYAAABTWwL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAADBYklEQVR4nOzdd9xcVYH/8e9JAgQwoZcAxkikSBEUFBDUIFUsuOqCrrr2surquupKrAFxwbY/dV0V1wL2vhaKQRAENHQp0pQWBCI9EAgBEs7vj5n75D43d2buzNxT7rmf9+vFi+e5U+6Zycz3Oe2eY6y1AgAAAAAAAFyYEroAAAAAAAAASBedTwAAAAAAAHCGzicAAAAAAAA4Q+cTAAAAAAAAnKHzCQAAAAAAAM7Q+QQAAAAAAABn6HwCAACIjDHmq8aYj4YuBwBgPOQ50GGstaHLAAAAgB6MMfMkfddau02f++wv6WOSniHpPmvtHC+FAwBURp6jzZj5BAAA0HwPSfqmpA+ELggAYCzkOZJE5xMAAIADxhhrjHlK7vcTjTHHdn+eZ4y51RjzPmPMncaYJcaYNxTva4xZX9JpkrYyxjzY/W+r4rmstRdaa78j6UYPLw0AWoU8B8ZH5xMAAEAYW0raQNLWkt4k6X+MMRvl72CtfUjSCyTdbq19Qve/2/0XFQDQB3kODEDnEwAAQBiPSTrGWvuYtfZUSQ9K2iFwmQAAwyPPgQHofAIAAAjjHmvtytzvyyU9IVRhAAAjI8+BAeh8AgAAcGO5pPVyv2854vOwNTEAhEWeA2Oi8wkAAMCNyyT9kzFmqjHmUEnPG/F57pC0iTFmg153MMZMMcZMl7RW51cz3Riz9ojnAwBMdpnIc2AsdD4BAAC48R5JL5a0VNKrJf1ilCex1l4r6QeSbjTGLC3bHUnScyU9LOlUSbO7P58+yvkAAGsgz4ExGWuZ+QcAAAAAAAA3mPkEAAAAAAAAZ+h8AgAAAAAAgDN0PgEAAAAAAMAZOp8AAAAAAADgDJ1PAAAAAAAAcGZa6AL4tummm9o5c+aELgYAYEw333yzyHMAaD7yHADScMkll9xtrd2s7LbWdT7NmTNHF198cehiAADGtOeee5LnAJAA8hwA0mCMWdzrNi67AwAAAAAAgDN0PgEAAAAAAMAZOp8AAAAAAADgDJ1PAAAAAAAAcIbOJwAAAAAAADhD5xMAAAAAAACcofMJAAAAAAAAztD5BAAAAAAAAGfofAIAAAAAAIAzdD4BAAAAAADAGTqfAAAAAAAA4AydTwAAAAAAAHCGzicAAAAAAAA4Q+cTAAAAAAAAnKHzCQAAAAAAAM7Q+QQAAAAAAABn6HwCAAAAAACAM3Q+AQAAAAAAwBk6nwAAAAAAAOBM4zufjDGHGmOuM8Zcb4w5KnR5AAAAAAAAsFqjO5+MMVMl/Y+kF0jaSdKrjDE7hS0VAAAAAAAAMtNCF2BMz5J0vbX2RkkyxvxQ0uGSrg5aKkRl+aMrJUlGpvN/M/n27Pfi7WbidlP4ffUxAAAAAADQX9M7n7aW9Lfc77dK2itQWRCpPY89Q8sfXeX0HFU6rDq/T76jKblt0HOp1+2mc2yKMZp/2FP1ij22Ge9FAQAAAABQg6Z3PlVijHmrpLdK0uzZswOXBr697+Ad9NiqxyVJ1naOWdlJv2estYX79X6MXf2gSb/n71+8b/G5tMbtq8/R6zE9y9z9/YcX/k1X3rqUzicgkCNPWCRJ+tHb9im9fdcFCyVJVy44pPZzZ89d5OJcZebOP0WSdMNxL6z9uecc1Xnum49/4cTPmZuPd3c+n+fJv7biuXodb5pe/3Zl/77ZbXV/rlx+Tot2XbBQy1Z0ZmDPmN6pdi9bsVIzpk+b9L10mQvS6lzK/Oht+wzMKpDnEnk+6nnIc/KcPF+TscWWbIMYY/aRtMBae0j39/mSZK09rtdj9txzT3vxxRd7KiHg3+7HnK7Dd9tKRx++S+iiAE7tueeeiiXPyyoCdcpXZHpValxXdqTVlbxeXDdSQhqmHP3uO87rieW9GFXo8ucbKb1+diVrtGQNFsl9J0I+l65e8oB2mjVz4vfYGijkOXnuE3k+vtDlb3OeZ7Jcjy3PjTGXWGv3LL2t4Z1P0yT9RdIBkm6TdJGkf7LWXtXrMXQ+IVWPP271+TP+oi/+7vqJY3/95Au01tRG7ysA9OSrsTJKQyQ/AuV6NKpfJWjUBozPkcQydY4W11VBDlHRLhtVLo6s13mefiP1PviYATF3/ila1a36TjV+P+M+OhSKyjqgYmuoSOR5hjwf7bnqKpNL5Dl5Pq6m5HmynU+SZIw5TNLnJU2V9E1r7Sf73Z/OJ6TqiK8u0oU33zvp2LEv3UWv2ftJgUoEuBXTSDkAYHTkOQCkoV/nU+PXfLLWnirp1NDlAEJ7w75zJnU+TZti9Oq9WOMM8Mn3NffZyGJ+lFxyf6nGKtsZZcyEGlF3pTha7Hv0uMo6Gq7Ol3/u0JdVhFL3iLbvEfJiDrm+jCxV5HkayPPJx8nzuJ5vkNTyvPEzn4bFzCekjjWf0BaMlANAGshzAEhD0jOfAAAYxzij2z7XAhlF3SN0+QVqs9Fyl6PkbR2pbbJB67uU3eZLce0b12vhDPr+Vf1+xpgtsSLPqyPPMQh5vhp5Xg9mPgGJuPehR3XYF87V3x9YIUk65wP7a/Ym6wUuFeAOI+UAkAbyHADS0G/mE9tgAYn44UW3THQ8SdKvr7g9YGmA9jjyhEWlW+Cmau78UwZu0V2HOUedssZaGW3h47Wn+v5W+Xz6+gyjmdqU6eS5e+T56Mjz9DDzCUiEtVaLbrhH//T1CzRvh8104hueFbpIgFOMlPsR6zbdobgsT2yvNRYhP4PDXOqU3TeTPSbEltxNQ577QZ6vyfXC3zG91hiQ5+nrN/OJzicgMSw4jrYI2VgZtNtIyGv6XVaMQjdcUhZ6DY3UG0q+1wcJ6cgTFunqJQ9op1kzG7OuCHlejjxvJvLcLfI8bnQ+5dD5hJStetxq7odO1a5bb6DPHbGbtt9iRugiAc64aKwM28ho08KRZZU71xW+XltUAxhvgezYsos894s8B+KSUp7T+ZRD5xNSVvzD/tO376M952wcqDSAWzFcpuFzRIqp3v3VNdJbVwOpankG7SYUqoGW4qUoxXVBXI6Sh/y+FhsisTVMypDnyCPP60Wej4c8Hw6dTzl0PiFlR3x1kS68+d6J369YcLBmTl8rYIkAd2JorMCNmCrqvspStlhsqjMD2jADothYKWu81NWgaUJjZBDyPF3kubycMxTyvPw+o0ohz+l8yqHzCaljzSe0RQyNlXGmSdetV8WnzhG7ufNP0SorTTWdUUYf6yyEqNiWnTN0Awpxqfq9Ki5am/E5gh46m6ogzycjz92ekzxHHnleLzqfcuh8QurofEJbxNBYGcRFJcHliFsV/RoodTde+o0eN7HxUFbmmEeNm/gelylrZOd/T0Fx0exMzA2UPPKcPM//3gTkeRjkefzofMqh8wmpuvGuB/X8z/1+0rHLPnaQNlxv7UAlAtxqQmMFADAYeQ4AaejX+TTNd2EAuPEPX/7jGsda1rcMBJFfpDbjcnSqytoDruVHUH1vcdxrZDmFUd0YLknxdd6yMrg+b36B2tAj5Sw4HSffeS5N/iyQ5+R5nef0dd6yMvg4bzbrSSLPm4KZT0AiVq56XJ/6zbX633Nv0vprT9Wfjz5ExpjQxQKcYaQcANJAngNAGpj5BLTAtKlT9OEX7qSfXHKrDt9tKzqeAA+afl1+P75HwGMQw6hxG6Qwq6GXcWcmjrOuUMp55EPK7x953pFi5oRGnvdGnq+JmU9AQpY/ulI7fawTjDcddxgdUEgaI+UAkAbyHADSwMwnoCWyjidJevlX/qifv2PfgKUB4BJrDLiRH8VNeUQ302/HpuLrruv9CLkjlO8ZIL2+pz6+v03Ykhsd5Lkb5Dl5XifyfHzMfAIScsoVS/TO718qSfr9B+bpSZusH7hEgDuxjZTnF6r1VTmIYfFxl2LZtjqmRku/soxazlCvL6b3tW6hF5EuU7yMI6ZGDHlOnvsuRwy
"text/plain": [
"<Figure size 1440x720 with 3 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"unit_id = 1\n",
"\n",
"fig, ax = plt.subplots(1,3, figsize=(20,10))\n",
"ax[0].plot(np.diff(stim['stim_ontime'], n=2), np.arange(20000-2))\n",
"ax[0].invert_yaxis()\n",
"ax[0].set_xlabel('$\\Delta^2[t_{on}]=\\Delta$ trial duration / s', size=20)\n",
"ax[0].set(ylabel='trials')\n",
"ax[1].eventplot(subseq_spike_times_locked['ontime'][unit_id], lineoffsets=1, linelengths=0.8)\n",
"ax[1].set(title='unit %d'%unit_id, yticks=[], xlabel='t/s (locked to trial onsets)')\n",
"ax[1].invert_yaxis()\n",
"ax[2].eventplot(subseq_spike_times_locked['offtime'][unit_id], lineoffsets=1, linelengths=0.8)\n",
"ax[2].set(title='unit %d'%unit_id, yticks=[], xlabel='t/s (locked to trial ends)')\n",
"ax[2].invert_yaxis()\n",
"plt.subplots_adjust(wspace=0)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"trial_spike_times = np.ndarray((num_unit, num_trial), dtype=object)\n",
"delta_t_on = 0.2 # 200ms\n",
"delta_t_off = 0.3\n",
"for trial_id in range(num_trial):\n",
" t_on, t_off = stim['stim_ontime'][trial_id]-delta_t_on, stim['stim_offtime'][trial_id]+delta_t_off\n",
" for unit_id in range(num_unit):\n",
" trial_spike_times[unit_id, trial_id] = spike_times[unit_id][np.where((spike_times[unit_id] < t_off) & (spike_times[unit_id] > t_on))] - stim['stim_ontime'][trial_id]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"dt = 0.001 # 1ms\n",
"T = delta_t_on + delta_t_off + trial_length\n",
"time_range = np.arange(0, T, dt) - delta_t_on\n",
"trial_onsets = time_range[np.where(np.diff(time_range//trial_length))]\n",
"num_bins = len(time_range)\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"PSTH = {}\n",
"\n",
"for key in key_list:\n",
" PSTH[key] = np.zeros((num_unit, num_stim[key], num_bins))\n",
" for stim_id, trials in enumerate(stim_id_trial[key]):\n",
" for unit_id in range(num_unit):\n",
" cumulative_spike_times = np.concatenate(trial_spike_times[unit_id, trials])\n",
" for t in cumulative_spike_times:\n",
" PSTH[key][unit_id, stim_id, int((t+delta_t_on)//dt)] += 1\n",
" # firing_rate[key][unit_id, stim_id] = np.histogram(cumulative_spike_times, dt_edges)\n",
" \n",
" PSTH[key] /= dt*len(trials)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"from scipy import signal\n",
"fs = 1/dt\n",
"b, a = signal.butter(5, 50, fs=fs, btype='low')\n",
"PSTH_filt = {}\n",
"for key in key_list:\n",
" PSTH_filt[key] = signal.filtfilt(b, a, PSTH[key], axis=-1)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import gridspec\n",
"\n",
"def plot_PSTH_raster(unit_id, key, stim_list=[], raster=True, plot=False, savepath=''):\n",
"\n",
" if len(stim_list)==0:\n",
" stim_iter = range(num_stim[key])\n",
" \n",
" else:\n",
" stim_iter = stim_list\n",
"\n",
" n_row = int(len(stim_iter)//4)+1\n",
" fig = plt.figure(figsize=(4*5, n_row*5))\n",
" gs = gridspec.GridSpec(n_row, 4)\n",
" for i, stim_id in enumerate(stim_iter):\n",
"\n",
" r, c = i //4, i%4\n",
"\n",
" if raster:\n",
" gs0 = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=gs[r, c], hspace=0)\n",
" ax0 = fig.add_subplot(gs0[0, :])\n",
"\n",
" ax0.eventplot(trial_spike_times[unit_id, stim_id_trial[key][stim_id]]+delta_t_on)\n",
" ax0.set(ylabel='trials', title=key_symbol[key] + ' = %d'%stim_val[key][stim_id])\n",
"\n",
" ax0.invert_yaxis()\n",
"\n",
"\n",
" ax1 = fig.add_subplot(gs0[1, :])\n",
" ax1.axvspan(0, trial_length, color='g', alpha=0.1)\n",
"\n",
" for t_on in trial_onsets:\n",
" ax1.axvline(t_on, c='k', alpha=0.5)\n",
"\n",
" ax1.plot(time_range, PSTH[key][unit_id, stim_id], alpha=0.3)\n",
" ax1.plot(time_range, PSTH_filt[key][unit_id, stim_id], 'r', label='smoothed PSTH')\n",
" ax1.set(xlabel='t / s')\n",
"\n",
" else:\n",
" ax1 = fig.add_subplot(gs[r, c])\n",
" ax1.axvspan(0, trial_length, color='g', alpha=0.1)\n",
" for t_on in trial_onsets:\n",
" ax1.axvline(t_on, c='k', alpha=0.5)\n",
" \n",
" ax1.plot(time_range, PSTH[key][unit_id, stim_id], alpha=0.3)\n",
" ax1.plot(time_range, PSTH_filt[key][unit_id, stim_id], 'r')\n",
" ax1.set(xlabel='t / s', title=key_symbol[key] + ' = %d'%stim_val[key][stim_id])\n",
" \n",
"\n",
" fig.suptitle('PSTH of unit %d' % unit_id, y=1.01)\n",
" plt.tight_layout()\n",
" \n",
" if plot:\n",
" plt.show()\n",
" \n",
" if len(savepath)>0:\n",
" plt.savefig(savepath, transparent=False, facecolor='white', bbox_inches='tight')\n",
" plt.close()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1\r"
]
}
],
"source": [
"for unit_id in range(1,2):\n",
" print(unit_id, end='\\r')\n",
" plot_PSTH_raster(unit_id=unit_id, key='orientation', stim_list=[], raster=True, plot=False, savepath='PSTH_orientation_unit%d'%unit_id)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"autocorr = np.convolve(B_spike[1].A[0][:100000], B_spike[1].A[0][::-1][:100000], mode='full')"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f882a32baf0>]"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAktUlEQVR4nO3deXxU1f3/8dchbAIKCCkii8F9XxEVxbqgdam1i7+6tErVfq1tbW2trfFr3Vr3WrXWflXUqq0bWnFrBJV9UYGwLyELECCQlUAWsifn98fchEkyk0ySmbkzd97PxyMP5p5779xP7kw+nHvuuecYay0iIhL/erkdgIiIhIcSuoiIRyihi4h4hBK6iIhHKKGLiHhE72gebPjw4TYlJSWahxQRiXvLly8vsdYmd7ZdVBN6SkoK6enp0TykiEjcM8ZsDWU7NbmIiHiEErqIiEcooYuIeIQSuoiIRyihi4h4hBK6iIhHKKGLiHiEEroIMHNdPiWVtW6HIdIjSuiS8Mqq67n1jRXc+Ooyt0MR6REldEl4jU2+SV7ydle5HIlIzyihi4h4hBK6iIhHKKGLiHiEErqIiEcooYuIeIQSuoiIRyihS0L675qdfLhyh9thiIRVVGcsEokVt721EoDvnjLK5UhEwkc1dPGsovIa0tbkh7y9jWAsItGghC6e9eNXlvDLt1ZQVdfQ4XbG+dcqo0ucU0IXz8rbXQ10nKittRjTuiynqJK9te3/E9i2q4o9VXXhDFEkrJTQJaG9s2x7u7LJT80POFDXuX+Zy+Sn5kcjLJFuUUKXhJZdWBmwfGluKQC1DY2tyksqVUOX2KWELhLE799bzVF/nElOUQX1jU1uhyPSKSV0kSDeW54HQEa+ErrEByV0ERGPUEIXEfGIkBK6Mea3xpj1xph1xpi3jTH9jTHjjDFLjDE5xphpxpi+kQ5WJByaZygCWLFtd6t1q7fvabf93dPXsqvNzdD/rtnJuwF6yIi4qdOEbowZBfwaGG+tPR5IAq4BHgeettYeDuwGbo5koCLhssovafu/ttZy5T8Wt9u+sraBRz7NaFV221sr+cP7ayIVoki3hNrk0hvYzxjTGxgA5AMXAP9x1r8OfDfs0YlEgXGeFdWDohLvOk3o1todwJPANnyJvAxYDuyx1jY/TpcHBBzlyBhzizEm3RiTXlxcHJ6oRSKgoqbjIQJEYl0oTS5DgSuBccDBwEDgklAPYK2daq0db60dn5yc3O1ARcJlbd6eVsvFlTXuBCISZqE0uUwGtlhri6219cB04GxgiNMEAzAa0ODSEvNyiip44JMNrcomP7WgS++Rt7sqnCGJhE0oCX0bcKYxZoAxxgAXAhuAucBVzjZTgI8iE6JI+LTtrRKq2oZ9Dxbtqapvt76xyVJT39iuXCSaQmlDX4Lv5ucKYK2zz1TgLuAOY0wOMAx4JYJxinRZOIfDnbOxqMP1d763mqPvnRm+A4p0Q0gzFllr7wfub1O8GZgQ9ohEwqzt8LiR8IGms5MYoCnoRHqguKKW5VtL3Q5DBFBCF+mR9NxSnvw80+0wRACN5SIeZgM8KhTuvuabigOPpy7iBiV08TzDvkb0VxZtCet7P/l5VljfT6Qn1OQinlNeU89HK3cE7OXy1eZd0Q9IJEqU0MVz7vlgHZ+s3ul2GCJRpyYX8ZzdewM/PPTO0m1hef9HZ2R0vpGIC5TQJSHsqaojdfrasLzX4hw120hsUpOLeEJ1XSN9kgxJvQyVta17slTU1gd8XF/Ea5TQxROOuW8m5x2VzMTDhrWatAJgwsOz3QlKJMqU0MUz5mUWU1Yd/Zq4icbYAiIhUBu6SA/lFLV/uGh7aRUpqWntrhZEIkkJXSQC5mf5Zud6N10TSUv0qMlFPGW1yzXilxduZnVeGaeOHeJqHJKYlNDFU5pcnun5oTRfH/V6vwkxRKJFTS4iIh6hhC4i4hFK6BK3Sipr2VhQ7nYYHSooq2k1xG52YQVF5TUuRiRepoQuceuCJ+dxyTML3Q6jQ3M2FnHhX+e3LF/09AImPKIHnSQylNAlbpWHebKKcGpw++6sJCQldJEImJVR6HYIkoCU0EVEPEIJXeJKSWUtKalpzM0sail7dna2ixGF5tj7ZpKSmtay/MMXv3IxGvEqJXSJK2t3lAHw2uLclrKnvoj9eT2r6hpbLS/dUupSJOJlSugSEwrKarjx1aVU1GjccpHuUkKXmPDsnGzmZhbz0SrNBSrSXUroEjdyiir5YMUOAFfGPQ+3qroGnp+3iUZ1cZQw0eBcElM6Sm0XPT0f62zghXHGn5iZyWtf5nLwkP5cefIot8MRD1ANXaJq555qsgor2pWHMueP9VhFtnnu04z8CraXVpGRX05hCMMCWGtZkFVMk2r20oZq6BJVEx+bA0DuY5e7HEnseGH+Jl6YvwmApF6GTY9c1uH2szKK+J9/pXP3pUfzs28eFo0QJU6ohi6u2FNVF3hFkGp4c23WS5oC/K5t29Obmiyle/edqz1VdeTtrgJgW2lVZAOUuKOELq44+U9ftFrubJ7l4+//LILRuGO6c4O3I3/9IpNT//wFJZW1gO+8PfjJhkiHJnFKCV1ignFa0dUq3Npn631jwuzeG+SKRsSPErrEhOYaun8rREVNPSmpaXywMs+doFxUU99ISmoaOUWVnW8s4lBCl5iVt7sagBfnb3Y5kuhrbmIR6YqQerkYY4YALwPH47sqvgnIBKYBKUAu8ENr7e5IBCmxb+a6fOZnFfPo90/s0n6vLt5CefW+G57WqaL//r3VvLfcVzNvTuyJbH5WMc/NzXE7DIlxodbQ/wbMtNYeDZwEZACpwGxr7RHAbGdZEtStb6zg7aXbu7zfg59s4OlZWe36oTcnc/BmD5eueigtQ8MiSKc6raEbYwYD5wI/AbDW1gF1xpgrgfOczV4H5gF3RSJIiS1zM4vol9SLiYcPD7j+o1U7OCx5EMePGhzye1bX+0YjnJaex942IxOKSGhCaXIZBxQDrxpjTgKWA7cDI6y1+c42BcCIQDsbY24BbgEYO3ZsjwMW99346jIg+MNBt7+zqsP1gbyb7quRZ+SXk5Ef2xM/i8SqUJpcegOnAs9ba08B9tKmecX6Gj4D9jiz1k611o631o5PTk7uabwinldYXsM6Z9x3ka4IpYaeB+RZa5c4y//Bl9ALjTEjrbX5xpiRQFHQdxCRkE16Yi51DU1uhyFxqNMaurW2ANhujDnKKboQ2AB8DExxyqYAH0UkQolZ1XWNFHUwmFRNfSNlVfuGua1tUNt4KMKZzJuaLEUV+z6jsup6qnWPwrNC7eXyK+BNY8wa4GTgEeAx4CJjTDYw2VmWBHLMfTOZ8MhsNhXve/hlxbZ9PVevfG4xJ/3p85bln76eHtX4xDc934SHZ7eM4njSg58z+an5LkclkRJSP3Rr7SpgfIBVF4Y1GolL/v3Es/2Gxs1sM0zuwuySqMUkPnM2+lpCiytqGXFAfwB27FG/fq/Sk6LSY1P+ubTTbWrq21/mZxa0Hxddoq+kspaU1DTS1uR3vrHENCV0iYo9Ve2njJuXqfvosSDL+Y/1ja+3uhyJ9JQmuJCwuuv9te3K3vh6K0u3lLYrV39zkfBSQpeI++OH6wKWf6hH2UXCSgldJE69uWQbyfv3Y92OcqZMPIRJR+jBvUSnhC4Sx56ZlQ3ArIxCzdMquikqnZufVcyGneVsLFCbdyxbsnlXy/yjC7N9n1mG85l1NJSAZonyDtXQpUObiitbdUtULTB2XT31a44cMYjXb5rA9a+07kqaOn0t10zoeHC8zuZ1ldinGrp0aK/GIo8rWYWVbMwP3L9/Z5AHiqyq6J6hhC7iMTe+tixg+cTH5nS4n2ro8U8JXTpk2s0lJCKxSm3oEpC1lnF3f8rkY77Rqtx/IC6JXxMfnU3K8IG89T9nuh2KhJESugTU3K46K6P14/kLsopdiEbCbWdZDTvLfCMwWvVz8Qw1uSSgytoGvv33hS3dEG98dSkz1+WzqbiSSU/M4ZzH5wSdYf7BTzZEM1SJorbNa+m5pXz//xZrso04ooSegL7MKWHdjnKe/Cw
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(autocorr)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def point_fft(pt, tapers, nfft):\n",
" \"\"\"\n",
" Compute Fourier spectra of a set of point processes given a set of tapers.\n",
" \n",
" Parameters\n",
" ----------\n",
" pt : ndarray\n",
" 2D array of point processes, columns should correspond to processes\n",
" tapers : ndarray\n",
" 2D array of tapers, rows should correspond to different tapers\n",
" nfft : int\n",
" number of points to include in the FFT, if nfft is greater than the number of time points the signals\n",
" will be zero-padded\n",
" \n",
" Returns\n",
" -------\n",
" out : ndarray\n",
" Fourier series of each tapered point process\n",
" \"\"\"\n",
" # allow 1D input\n",
" if pt.ndim == 1:\n",
" pt = pt[:, np.newaxis]\n",
" \n",
" assert pt.shape[0] == tapers.shape[1], \"Tapers and point-processes must have same number of time points.\"\n",
" \n",
" # apply tapers to point process\n",
" tapered_pt = np.array([pt.T * taper for taper in tapers])\n",
" \n",
" # compute mean event rate (repeated for each taper)\n",
" mean_rates = np.ones((tapers.shape[0], pt.shape[-1])) * pt.mean(axis=0)\n",
" \n",
" # compute FT of tapers\n",
" taper_fft = fft.fft(tapers, nfft)\n",
" \n",
" # compute a 'baseline' spectrum by scaling the taper FT by the mean rate\n",
" baseline = np.array([mean_rates.T * taper for taper in taper_fft.T]).T\n",
" \n",
" # compute FT of tapered point process and subtract baseline \n",
" return fft.fft(tapered_pt, nfft) - baseline "
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"from scipy import io, fft, signal\n",
"from scipy.signal import windows, detrend\n",
"\n",
"def mtptcoh(pt, ct, fs=1, nperseg=None, nfft=None, noverlap=None, nw=3, ntapers=None, detrend_method='constant',\n",
" nspike_min=4, progress_report=True):\n",
" \"\"\"\n",
" Compute the coherence between a set of point processes and a set of continuous signals. \n",
" \n",
" Parameters\n",
" ----------\n",
" pt : ndarray\n",
" 2D array of point processed, columns should correspond to processes (termed units)\n",
" ct : ndarray\n",
" 2D array of continous signals, columns should correspond to channels\n",
" fs : float (default = 1)\n",
" sampling frequency\n",
" nperseg : int, None (default = None)\n",
" number of data points per segment, if None nperseg is set to 256\n",
" nfft : int, None (default = None)\n",
" number of points to include in scipy.fft.fft, if None nfft is set to 2 * nperseg, if nfft > nperseg data \n",
" will be zero-padded\n",
" noverlap : int, None (default = None)\n",
" amout of overlap between consecutive segments, if None noverlap is set to nperseg / 2\n",
" nw : int (default = 3)\n",
" time-frequency bandwidth for Slepian tapers, passed on to scipy.signal.windows.dpss\n",
" ntapers : int, None (default = None)\n",
" number of tapers, passed on to scipy.signal.windows.dpss, if None ntapers is set to nw * 2 - 1 (as \n",
" suggested by original authors)\n",
" detrend_method : {'constant', 'linear'} (default = 'constant')\n",
" method used by scipy.signal.detrend to detrend each segment\n",
" nspike_min : int\n",
" smallest number of discrete events (per time segment) allowed for a unit to be considered in the\n",
" computation\n",
" \n",
" Returns\n",
" -------\n",
" f : ndarray\n",
" frequency bins\n",
" Sxx_pt, Sxx_ct : ndarray x 2\n",
" 'auto-spectra' of point processes and continuos signals respectively\n",
" Sxy : ndarray\n",
" cross-spectra between every combination of point process and continuous signal channels\n",
" Cxy : ndarray\n",
" coherence magnitude between every combination of point process and continuous signal channels\n",
" phloc : ndarray\n",
" phase locking between every combination of point process and continuous signal channels, phase locking\n",
" index in each frequency bin can be computed by taking the angle of matrix each entry\n",
" \"\"\"\n",
" # allow 1D inputs\n",
" if pt.ndim == 1:\n",
" pt = pt[:, np.newaxis]\n",
" \n",
" if ct.ndim == 1:\n",
" ct = ct[:, np.newaxis]\n",
" \n",
" assert pt.shape[0] == ct.shape[0], \"Signals must be same length.\"\n",
" \n",
" # set some default for parameters values\n",
" if nperseg is None:\n",
" nperseg = 256\n",
" \n",
" if nfft is None:\n",
" nfft = nperseg * 2\n",
" \n",
" if noverlap is None:\n",
" noverlap = nperseg / 2\n",
" \n",
" if ntapers is None:\n",
" ntapers = 2 * nw - 1\n",
" \n",
" stepsize = nperseg - noverlap\n",
" nsegs = int(np.floor(pt.shape[0] / stepsize))\n",
"\n",
" fftnorm = np.sqrt(2 / nfft) # value taken from original matlab function\n",
" \n",
" # initialize auto-spectra and cross-spectra arrays\n",
" Sxx_pt = np.zeros((pt.shape[1], nfft), dtype='complex128')\n",
" Sxx_ct = np.zeros((ct.shape[1], nfft), dtype='complex128')\n",
" Sxy = np.zeros((pt.shape[1], ct.shape[1], nfft), dtype='complex128')\n",
" phloc = np.zeros((pt.shape[1], ct.shape[1], nfft), dtype='complex128')\n",
" \n",
" # get FFT frequency bins\n",
" f = fft.fftfreq(nfft, 1 / fs)\n",
" \n",
" # get tapers\n",
" tapers = windows.dpss(nperseg, nw, Kmax=ntapers)\n",
"\n",
" # keep track of number of segments used for each unit \n",
" unit_nsegs = np.zeros(pt.shape[1]) \n",
" # loop over segments\n",
" for seg_ind in range(nsegs):\n",
" if progress_report:\n",
" print('%d'%(100*(seg_ind+1)/nsegs)+'%', end='\\r')\n",
" \n",
" # prepare segment\n",
" i0 = int(seg_ind * stepsize)\n",
" i1 = int(seg_ind * stepsize + nperseg)\n",
" if i1 > len(pt): # stop if segment is out of range of data\n",
" break\n",
" seg_pt = pt[i0:i1, :]\n",
" seg_ct = ct[i0:i1, :]\n",
" \n",
" # get nspikes\n",
" nspikes = seg_pt.sum(axis=0)\n",
" valid_units = np.where(nspikes >= nspike_min)[0]\n",
" \n",
" # restrict analysis for this segment to units with sufficient spikes\n",
" seg_pt = seg_pt[:, valid_units]\n",
" unit_nsegs[valid_units] += 1\n",
" \n",
" # compute FT of point process data\n",
" pxx_pt = point_fft(seg_pt, tapers, nfft) / fftnorm\n",
" \n",
" # compute 'auto-spectra' of point process\n",
" Sxx_pt[valid_units] += (pxx_pt * np.conjugate(pxx_pt)).sum(axis=0) \n",
" \n",
" # detrend continuous data\n",
" seg_ct = detrend(seg_ct, type=detrend_method, axis=0)\n",
"\n",
" # apply tapers\n",
" tapered_seg_ct = np.array([seg_ct.T * taper for taper in tapers]) \n",
"\n",
" # compute FFT for each channel-taper combination\n",
" pxx_ct = fft.fft(tapered_seg_ct, nfft) / fftnorm\n",
" \n",
" # compute 'auto-spectra' of continuous process\n",
" Sxx_ct += (pxx_ct * np.conjugate(pxx_ct)).sum(axis=0)\n",
" \n",
" # for each unit-channel combination\n",
" for unit in range(seg_pt.shape[1]):\n",
" for channel in range(seg_ct.shape[1]):\n",
" # compute cross-spectra, averaging over tapers\n",
" Sxy[unit, channel] += (pxx_pt[:, unit, :] * np.conjugate(pxx_ct[:, channel, :])).sum(axis=0)\n",
" # take only phase values and compute PLI, averagins over tapers\n",
" pt_phase = np.angle(pxx_pt[:, unit, :])\n",
" ct_phase = np.angle(pxx_ct[:, channel, :])\n",
" phloc[unit, channel] += (np.exp(1j * pt_phase) * np.conjugate(np.exp(1j * ct_phase))).sum(axis=0) \n",
" \n",
" # get normalization factors to account for summing over tapers and segments\n",
" unit_csdnorm = unit_nsegs * ntapers\n",
" continuous_csdnorm = nsegs * ntapers\n",
" \n",
" # apply normalization factors (divide by zero will result in NaN)\n",
" Sxx_pt = (Sxx_pt.T / unit_csdnorm).T\n",
" Sxx_ct = (Sxx_ct.T / continuous_csdnorm).T\n",
" Sxy = (Sxy.T / unit_csdnorm).T\n",
" phloc = (phloc.T / unit_csdnorm).T\n",
" \n",
" # compute power normalization matrix\n",
" powernorm = np.zeros((pt.shape[1], ct.shape[1], nfft))\n",
" for unit in range(pt.shape[1]):\n",
" for channel in range(ct.shape[1]):\n",
" powernorm[unit, channel] = np.sqrt((np.abs(Sxx_pt[unit])) * (np.abs(Sxx_ct[channel]))) \n",
" \n",
" # apply normalization to cross-spectra to get coherence\n",
" Cxy = np.abs(Sxy) / powernorm\n",
" \n",
" return f, Sxx_pt, Sxx_ct, Sxy, Cxy, phloc\n",
"\n",
"def inds2train(inds_list, n):\n",
" \"\"\"\n",
" Convert a set of event indices (time points) into a binary time-series. It is assumed that all channels in\n",
" the input set occur over the same period and that the indices refer to the same time base.\n",
" \n",
" Parameters\n",
" ----------\n",
" inds_list : list, dict, ndarray\n",
" sequence containing indices (time points) at which events occur for each of a set of \n",
" point processes\n",
" n : int\n",
" total length of segment in which events occur\n",
" \n",
" Returns\n",
" -------\n",
" event_train : ndarray\n",
" binary time-series array for each member of input set, columns correspond to channels\n",
" \"\"\"\n",
" event_train = np.zeros((n, len(inds_list)), dtype='uint8')\n",
" for ch in range(len(inds_list)): \n",
" event_train[inds_list[ch], ch] = 1\n",
" \n",
" return event_train\n",
"\n",
"def point_fft(pt, tapers, nfft):\n",
" \"\"\"\n",
" Compute Fourier spectra of a set of point processes given a set of tapers.\n",
" \n",
" Parameters\n",
" ----------\n",
" pt : ndarray\n",
" 2D array of point processes, columns should correspond to processes\n",
" tapers : ndarray\n",
" 2D array of tapers, rows should correspond to different tapers\n",
" nfft : int\n",
" number of points to include in the FFT, if nfft is greater than the number of time points the signals\n",
" will be zero-padded\n",
" \n",
" Returns\n",
" -------\n",
" out : ndarray\n",
" Fourier series of each tapered point process\n",
" \"\"\"\n",
" # allow 1D input\n",
" if pt.ndim == 1:\n",
" pt = pt[:, np.newaxis]\n",
" \n",
" assert pt.shape[0] == tapers.shape[1], \"Tapers and point-processes must have same number of time points.\"\n",
" \n",
" # apply tapers to point process\n",
" tapered_pt = np.array([pt.T * taper for taper in tapers])\n",
" \n",
" # compute mean event rate (repeated for each taper)\n",
" mean_rates = np.ones((tapers.shape[0], pt.shape[-1])) * pt.mean(axis=0)\n",
" \n",
" # compute FT of tapers\n",
" taper_fft = fft.fft(tapers, nfft)\n",
" \n",
" # compute a 'baseline' spectrum by scaling the taper FT by the mean rate\n",
" baseline = np.array([mean_rates.T * taper for taper in taper_fft.T]).T\n",
" \n",
" # compute FT of tapered point process and subtract baseline \n",
" return fft.fft(tapered_pt, nfft) - baseline "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spike triggered covariance\n",
"How can we define spike triggered average/covariance for a stimulus with parameters on a circle/[klein bottle]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.5 64-bit",
"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.10.5"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}