Mikkel
/
septum-mec
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
septum-mec/actions/stimulus-spike-lfp-response.../data/10-calculate-stimulus-spike...

893 lines
215 KiB
Plaintext
Raw Normal View History

actions/stimulus-spike-lfp-response-other-drive/
2021-03-10 12:30:22 +00:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculate spike-lfp coherence using other drive"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"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 scipy\n",
"import scipy.signal\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",
"import elephant as el\n",
"import neo\n",
"from scipy.signal import find_peaks\n",
"from scipy.interpolate import interp1d\n",
"from matplotlib import mlab\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": 6,
"metadata": {},
"outputs": [],
"source": [
"data_loader = dp.Data()\n",
"actions = data_loader.actions\n",
"project = data_loader.project"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"output = pathlib.Path('output/stimulus-spike-lfp-response-other-drive')\n",
"(output / 'data').mkdir(parents=True, exist_ok=True)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"identify_neurons = actions['identify-neurons']\n",
"# sessions = pd.read_csv(identify_neurons.data_path('sessions'))\n",
"units = pd.read_csv(identify_neurons.data_path('units'))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def get_lim(action_id):\n",
" stim_times = data_loader.stim_times(action_id)\n",
" if stim_times is None:\n",
" return [0, np.inf]\n",
" stim_times = np.array(stim_times)\n",
" return [stim_times.min(), stim_times.max()]\n",
"\n",
"def get_mask(lfp, lim):\n",
" return (lfp.times >= lim[0]) & (lfp.times <= lim[1])\n",
"\n",
"def zscore(a):\n",
" return (a - a.mean(0)) / a.std(0)\n",
"\n",
"def compute_stim_freq(action_id):\n",
" stim_times = data_loader.stim_times(action_id)\n",
" if stim_times is None:\n",
" return np.nan\n",
" stim_times = np.array(stim_times)\n",
" return 1 / np.mean(np.diff(stim_times))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def signaltonoise(a, axis=0, ddof=0):\n",
" a = np.asanyarray(a)\n",
" m = a.mean(axis)\n",
" sd = a.std(axis=axis, ddof=ddof)\n",
" return np.where(sd == 0, 0, m / sd)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def compute_energy(p, f, f1, f2):\n",
" if np.isnan(f1) or np.all(np.isnan(p)):\n",
" return np.nan\n",
" mask = (f > f1) & (f < f2)\n",
" df = f[1] - f[0]\n",
" return np.sum(p[mask]) * df"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def find_theta_peak(p, f, f1, f2):\n",
" if np.all(np.isnan(p)):\n",
" return np.nan, np.nan\n",
" mask = (f > f1) & (f < f2)\n",
" p_m = p[mask]\n",
" f_m = f[mask]\n",
" peaks, _ = find_peaks(p_m)\n",
" idx = np.argmax(p_m[peaks])\n",
" return f_m[peaks[idx]], p_m[peaks[idx]]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def compute_half_width(p, f, m_p, m_f):\n",
" if np.isnan(m_p):\n",
" return np.nan, np.nan\n",
" m_p_half = m_p / 2\n",
" half_p = p - m_p_half\n",
" idx_f = np.where(f <= m_f)[0].max()\n",
" idxs_p1, = np.where(np.diff(half_p[:idx_f + 1] > 0) == 1)\n",
" if len(idxs_p1) == 0:\n",
" return np.nan, np.nan\n",
" m1 = idxs_p1.max()\n",
" idxs_p2, = np.where(np.diff(half_p[idx_f:] > 0) == 1)\n",
" m2 = idxs_p2.min() + idx_f\n",
" assert p[m1] < m_p_half < p[m1+1], (p[m1], m_p_half, p[m1+1])\n",
" assert p[m2] > m_p_half > p[m2+1], (p[m2], m_p_half, p[m2+1])\n",
" \n",
" f1 = interp1d([half_p[m1], half_p[m1 + 1]], [f[m1], f[m1 + 1]])(0)\n",
" f2 = interp1d([half_p[m2], half_p[m2 + 1]], [f[m2], f[m2 + 1]])(0)\n",
" return f1, f2"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# p = np.load('debug_p.npy')\n",
"# f = np.load('debug_f.npy')\n",
"# compute_half_width(p, f, 0.01038941, 30.30187709636872)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# plt.plot(f, p)\n",
"# plt.xlim(29.9,30.6)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def compute_stim_peak(p, f, s_f):\n",
" if np.isnan(s_f):\n",
" return np.nan\n",
" return interp1d(f, p)(s_f)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def compute_spike_lfp_coherence(anas, sptr, NFFT):\n",
"\n",
" sigs, freqs = el.sta.spike_field_coherence(anas, sptr, **{'nperseg': NFFT})\n",
" return sigs, freqs"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def butter_bandpass(lowcut, highcut, fs, order=5):\n",
" nyq = 0.5 * fs\n",
" low = lowcut / nyq\n",
" high = highcut / nyq\n",
" b, a = scipy.signal.butter(order, [low, high], btype='band')\n",
" return b, a\n",
"\n",
"\n",
"def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n",
" b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n",
" y = scipy.signal.filtfilt(b, a, data)\n",
" return y\n",
"\n",
"# def compute_spike_phase_func(lfp, times, return_degrees=False):\n",
"# x_a = hilbert(lfp)\n",
"# x_phase = np.angle(x_a)\n",
"# if return_degrees:\n",
"# x_phase = x_phase * 180 / np.pi\n",
"# return interp1d(times, x_phase)\n",
"\n",
"\n",
"def vonmises_kde(data, kappa=100, n_bins=100):\n",
" from scipy.special import i0\n",
" bins = np.linspace(-np.pi, np.pi, n_bins)\n",
" x = np.linspace(-np.pi, np.pi, n_bins)\n",
" # integrate vonmises kernels\n",
" kde = np.exp(kappa * np.cos(x[:, None] - data[None, :])).sum(1) / (2 * np.pi * i0(kappa))\n",
" kde /= np.trapz(kde, x=bins)\n",
" return bins, kde\n",
"\n",
"\n",
"def spike_phase_score(phase_bins, density):\n",
" import math\n",
" import pycircstat as pc\n",
" ang = pc.mean(phase_bins, w=density)\n",
" vec_len = pc.resultant_vector_length(phase_bins, w=density)\n",
" # ci_lim = pc.mean_ci_limits(head_angle_bins, w=rate)\n",
" return ang, vec_len"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def compute_clean_lfp(anas, width=500, threshold=2):\n",
" anas = np.array(anas)\n",
" idxs, = np.where(abs(anas) > threshold)\n",
" for idx in idxs:\n",
" anas[idx-width:idx+width] = 0 # TODO AR model prediction\n",
" return anas, idxs\n",
"\n",
"\n",
"def compute_clean_spikes(spikes, idxs, times, width=500):\n",
"\n",
" for idx in idxs:\n",
" t0 = times[idx-width]\n",
" stop = idx + width\n",
" if stop > len(times) - 1:\n",
" stop = len(times) - 1 \n",
" t1 = times[stop]\n",
" mask = (spikes > t0) & (spikes < t1)\n",
" spikes = spikes[~mask]\n",
" spikes = spikes[spikes <= times[-1]]\n",
" return spikes\n",
"\n",
"\n",
"def prepare_spike_lfp(anas, sptr, t_start, t_stop):\n",
"\n",
" t_start = t_start * pq.s if t_start is not None else 0 * pq.s\n",
" sampling_rate = anas.sampling_rate\n",
" units = anas.units\n",
" times = anas.times\n",
" if t_start is not None and t_stop is not None:\n",
" t_stop = t_stop * pq.s\n",
" mask = (times > t_start) & (times < t_stop)\n",
" anas = np.array(anas)[mask,:]\n",
" times = times[mask]\n",
" \n",
" # take best channel from other drive\n",
" best_channel = np.argmax(signaltonoise(anas))\n",
"# best_channel = np.random.choice(anas.shape[1])\n",
" \n",
" cleaned_anas, idxs = compute_clean_lfp(anas[:, best_channel])\n",
" \n",
" cleaned_anas = neo.AnalogSignal(\n",
" signal=cleaned_anas * units, sampling_rate=sampling_rate, t_start=t_start\n",
" )\n",
" \n",
" spike_units = sptr.units\n",
" spike_times = sptr.times\n",
" spike_times = compute_clean_spikes(spike_times, idxs, times)\n",
"\n",
" sptr = neo.SpikeTrain(\n",
" spike_times[(spike_times > t_start) & (spike_times < times[-1])], units=spike_units,\n",
" t_start=t_start, t_stop=times[-1]\n",
" )\n",
"\n",
" return cleaned_anas, sptr, best_channel"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"def compute_spike_phase_func(lfp, times, return_degrees=False):\n",
" from scipy.fftpack import next_fast_len\n",
" x_a = scipy.signal.hilbert(\n",
" lfp, next_fast_len(len(lfp)))[:len(lfp)]\n",
"# x_a = hilbert(lfp)\n",
" x_phase = np.angle(x_a, deg=return_degrees)\n",
" return interp1d(times, x_phase)\n",
"\n",
"\n",
"def compute_spike_phase(lfp, spikes, flim=[6,10]):\n",
" \n",
" sample_rate = lfp.sampling_rate.magnitude\n",
" \n",
" # sometimes the position is recorded after LFP recording is ended\n",
" times = np.arange(lfp.shape[0]) / sample_rate\n",
" \n",
" spikes = np.array(spikes)\n",
" spikes = spikes[(spikes > times.min()) & (spikes < times.max())]\n",
" \n",
" filtered_lfp = butter_bandpass_filter(\n",
" lfp.magnitude.ravel(), *flim, fs=sample_rate, order=3)\n",
"\n",
" spike_phase_func = compute_spike_phase_func(filtered_lfp, times)\n",
" \n",
" return spike_phase_func(spikes), filtered_lfp"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, 100)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAIHCAYAAACBqzFTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd2Bb1fn/8Y9sZw9nJ2SHDIfsbbtQ4EsZpeXHKNDSlrJLGWWWFkqhmL1byt4jhDDCLgQIJUBI8Mh0thJnx44TO8OJ95B+f8iSJetKutqy/X79E/vq3qsTWffe85zznHMsdrtdAAAAAABEW1K8CwAAAAAAaBsIQAEAAAAAMUEACgAAAACICQJQAAAAAEBMEIACAAAAAGKCABQAAAAAEBMpsXqj9PR0+6BBg2L1dgAAAACAGOnZs6cWL178ldVq/bm//WIWgA4aNEgffvhhrN4OAAAAABBDaWlpfQLtQwouAAAAACAmCEABAAAAADFBAAoAAAAAiAkCUAAAAABATBCAAgAAAABiIuRZcNPS0v4u6UxJ7SU9a7VaX4lYqQAAAAAArU5IPaBpaWknSvqJpGMlnSBpSATLBAAAAABohULtAT1N0hpJH0nqLumvESsRAAAAAKBVCnUMaB9JMySdL+kqSW+lpaVZIlYqAAAAAECrE2oP6H5JG61Wa60ka1paWrWkvpL2RaxkAAAAAIBWJdQe0MWSfp6WlmZJS0sbKKmLHEEpAAAAAACGQgpArVbrZ5JWSsqT9F9J11qt1oZIFgwAAAAA0LqEvAyL1Wr9WyQLAgAAAABo3UJNwQUAAAAAICgEoAAAAACAmCAABQAAAADEBAEoAAAAACAmCEABAAAAADFBAAoAAAAAiAkCUACAoR37KzT8ts+1Yc/heBcFAAC0EgSgAABDX64tliR9tLIwziUBAACtBQEoAAAAACAmCEABAAAAADFBAAoAAAAAiAkCUACAX3a7Pd5FAAAArQQBKADAkMUS7xIAAIDWhgAUAAAAABATBKAAAAAAgJggAAUA+MUQUAAAECkEoAAAQxYxCBQAAEQWASgAAAAAICYIQAEAAAAAMUEACgAAAACICQJQAIBfzEEEAAAihQAUAGDIwhxEAAAgwghAAQAAAAAxQQAKAAAAAIgJAlAAgF92BoECAIAIIQAFAAAAAMQEASgAAAAAICYIQAEAAAAAMUEACgAA0IYU7DsiO4O7AcQJASgAwC+7qKgCrcXCjXt18r8W6ZNVRfEuCoA2igAUAGDIYrHEuwgAImzT3nJJ0oY9h+NcEgBtFQEoAAAAACAmCEABAIYqa+olSXsPV8e5JAAAoLUgAAUAGPps9R5J0vw1xXEuCQAAaC0IQAEAhhgCCgAAIo0AFABgiFUaAABApBGAAgAAAABiggAUAAAAABATBKAAAEMTB6dKkvp26xDnkgAAgNaCABQAYGjiIEcAevqEAXEuCQAAaC0IQAEAANoIJhcDEG8EoAAAAG0NyywBiBMCUAAAIEmqrbfpwfkbdKS6Lt5FAQC0UgSgAAC/SNlrO+Yt36UXFm3Vv7/eHO+iAABaKQJQAIAhCyl6bU59g6O1od5mi3NJAACtFQEoAAAAgDbDbrdrbWFZvIvRZhGAAgAAZX26Tnd9ui7exQCAqJuTu1NnPLVYizaVxLsobRIBKAAA0Os/bo93EQAgJjbuOSxJ2nGgMs4laZsIQAEAAAAAMUEACgAAAACICQJQAADgYXb2jngXAQDQShGAAgAAtBF2sbAvgPgiAAUAAGhjLGKhXwDxQQAKAPCLHhMAABApBKAAAEP0jwChq6pt0Ms/bJXNRgMOEEt7yqp0zVvLVVXbEO+iwAcCUAAAgAh7fIFV932+Qf9dXRTvogBtygPzN2r+mmItWF8c76LABwJQAACACCurqpMk1dTZ4lwSAEgsBKAAAAAA2h47KfLxQAAKADDEYxkA0BpZmOQgrghAAQB+sVwDAKAlW7CuWMVl1fEuBhoRgAIAAABota58c7l+9eySeBcDjQhAAQAAALQKvnJ2iugBTRgEoAAAAG0Ec64AiDcCUACAX3amIwJaHSZhARAvBKAAAEPUTwEAQKQRgAIAAAAAYoIAFAAAIMJIXAfii/HOiYsAFAAAIFrIZQdiivHNiY8AFAAAAECbQydpfBCAAgAAAGgzLKQmxBUBKAAAAAAgJsIKQNPS0vqlpaXtSktLGxupAgEAEsucnJ36ZFVhvIsBAABagZAD0LS0tHaSXpBUFbniAAAS0dMLC+JdBAAATLMzwjNhhdMD+pik5yUVRagsAAAAABAyM6M7CU7jK6QANC0t7RJJJVar9avIFgcAAADRxhQsANdBvITaA3qZpFPS0tK+kzRF0uy0tLQBESsVAAAAAKDVSQnlIKvVerzz58Yg9Cqr1VocqUIBAAC0ZHYy/ICEYOdiTDgswwIAMGYhOQkIlXOMGVcRAHgKqQfUndVqPTEC5QAAJDDaj4HQWGjIAeKCjs/ERQ8oAAAAgFYhmEYfYtT4IAAFAAAA0GZYSI6PKwJQAAAAAEBMEIACAAAAAGKCABQAAKCNYEkKAPFGAAoAANDGMDkvWjvaWhIXASgAICB6TQAALUHzthUeX4mHABQAAABAq0avf+IgAAUAAAAAxAQBKAAAQKSR9gcAhghAAQAAooSsPwDwRAAKAAAQQfd9tl4friyMdzGANs1MEgITFMUHASgAwBA9N0BoXl68Ld5FANouEw8vJiSKLwJQAEBANBIDAIBIIAAFAABoIV5ZvE3nPfdjvIsBtBg0oCaelHgXAAAAAObc+9n6sI5nzBvaKrJuEwc9oAAAAG2Mheo4gDghAAUAAADQqtgbu/vtdPsnHAJQAAAAAK2Cr959C1PfJgwCUAAAAABATBCAAgACI4MJCAmdLkBs2XlgJTwCUAAAAACtCim3iYsAFAAAIEqoAwPxRX9o4iEABQAAANCqNJ/9lragxEEACgAAAKBVYI3bxEcACgAwROpg28baeQCAaCAABQAAaCNoVgAQbwSgAAAALczO/ZVhHU+GA4B4IQAFAABoYUorauJdBAAICQEoAAAAgFaFdPPERQAKAAAAoFUgvTzxEYACAAwxCSoAoKXjWZZ4CEABAAAAtGpGPaMsNxUfBKAAgIB4RAOhsSg6+YDUm4HQkaUbXwSgAAAALQ4RKOCX3fkP10qiIQAFAAAA0Cr46t2MVjYCgkcACgAAAACICQJQAACANoKxowDijQAUAACgjSEZEUC8EIACAAwFs5h30aEqVdc1RK8wiDl6ygC0BtzLEg8BKAAgoEBrpf3koYW65q0VMSoNACrVgH9es9/S7Z8wCEABABGxcOO+eBcBaDOIPwFjwWTvID4IQAEAAKKEyjAAeCIABQAAANDmkEkQHwSgAAAAANoMC6kJcUUACgAAAACICQJQAACAFoZZcAH/uEYSFwEoACAgnuMAgJbAwnorCY8AFAAAoIUJtDavz+NoTkIbRViaOAhAAQAA2homYUEbFmoDDiKDABQAAABAm0MzTHwQgAIAAAAAYoIAFABgKJSJHA5X10WhJIgHEtQSG38fwD/nNUK2beIhAAUAhMV9LM2krAVxLAkAoK3zNbyZYc+JgwAUABAQLcgAgNaGR1t8EIACAABEiSVK3S40CgHmGC09FK3rEuYQgAIAAABo1UKZ1wDRQQAKAAgLPTEAAMAsAlAAAAAArQqNo4mLABQAEBae8UDsGY1rM3UcFyxaOYZ3Jj4CUAAAgDaGOjqAeCEABQAAAADEBAEoAABAS0MqLYAWigAUAAAgSkh1BQBPBKAAAAAAWiUm3ko8BKAAgLDYebq3SvxdExt/HSA4RrPjcpuLDwJQAIAhs1PZ8/wGfGNJCADwRAAKAAAAAIgJAlAAQEChLnoPAADgjgAUAACghWHsGoCWigAUAACgjSBuRVvhzNzhO594CEABAABamNqGhrCOZ3IktF7GX26+8okjJZSD0tLS2kl6VdJwSR0k3We1Wj+NYLkAAC0EqYBA7JVV1cW7CECLwJJSiSfUHtALJe23Wq0/lfRzSU9HrkgAAAA
"text/plain": [
"<Figure size 1152x648 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA74AAAIUCAYAAAA9jvA9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5wdZ3X4/8/M3LZ7tzetVr15Jau4YGMbm2IMBpMYEvwDEkyHAIEvJQRCSEiBEIhDSIAABtNb6KaDMTY27kWWLUuWtOraKm2/u3v7lN8fz8wtqivtap/11Xm/Xvu6c/vcnVvmzDnPeQzP8xBCCCGEEEIIISqVqXsFhBBCCCGEEEKIs0kCXyGEEEIIIYQQFU0CXyGEEEIIIYQQFU0CXyGEEEIIIYQQFU0CXyGEEEIIIYQQFU0CXyGEEEIIIYQQFS2kewVmy6WXXuotWbJE92oIMSO5XI5IJKJ7NYSYMXkvi0oh72VRCeR9LCrFU089NdzV1dV6JvetmMC3tbWVW2+9VfdqCDEjO3fuZN26dbpXQ4gZk/eyqBTyXhaVQN7HolJ0dnYeOtP7SqmzEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKJoGvEEIIIYQQQoiKds4Gvt0jKa79nz/ymq88rHtVhBBCCCGEEEKcRSHdK6DLV+7bz+4jU+w+MkXecQlb5+wxACGEEEIIIYSoaOdstHdgOFlY7h5NaVwTIYQQQgghhBBn0zkb+B4cSbK0qRqAfYNTmtdGCCGEEEIIIcTZorXUubOzsx74DlAHRID3dXV1PehfZwE/AL7S1dV122w+b8526RtL8/pnLefr9x9k31Dy1HcSQgghhBBCCPG0pDvj+z7gzq6urucCbwA+D9DZ2bkKuAe49Gw8afdoCteDjYvqaauNsm9IMr5CCCGEEEIIUal0N7f6HyDrL4eAjL9cA7wF+ODZeNJDIyrDu7wlzvKWeOG8EEIIIYQQQojKM2eBb2dn55uBvznq4jd2dXU92tnZ2Y4qeX4vQFdX11b/PtN+fM/z2Llz57Ru++jOBADZ4T4ibpZD47lp31eIsymTych7UVQEeS+LSiHvZVEJ5H0sxBwGvl1dXV8Fvnr05Z2dnRuB7wPv7+rq+uOZPr5hGKxbt25at72tdzcwwjMvXM8vD21n5/Dhad9XiLNp586d8l4UFUHey6JSyHtZVAJ5Hwuhv7nV+cCPgFcFWd65kEjnqY2FsEyDxuow4+k8ruthmsZcrYIQQgghhBBCiDmiu7nVJ4AY8JnOzs67Ozs7fz4XT5pI52moDgPQWB3BcT0mM3bhes/zGJzInOjuQgghhBBCCCGeRrRmfLu6ul52iuvfcDaeN5HOU19VDHwBxlI56v1g+KeP9/HBnzzJH/72eSzx5/oVQgghhBBCCPH0pDvjq8V4KlcMfOPqdCyVK1x/755h8o7H7TuOaFk/IYQQQgghhBCz55wMfBPpPA1VKtPb4Gd8x1P5wvWPHhwF4PanDs/9ygkhhBBCCCGEmFXnbOBb52d8m/zAdzSpMr6HExl6x9I0Vod59OAombyjbT2FEEIIIYQQQszcORf4ep53THMrKJY6b+9Tc/y+9IIOXA8OjaTK7t91eJIX/c899IyWXy6EEEIIIYQQYn6q+MA3mbVJ54pZ21TOIe94hTG+tbEQplEsdR7wuzk/57xWAA4MJ8se7wM/3krXkUnu2TM0F6svhBBCCCGEEGKGKjrw9TyPV93yIG/7zmOFyxJpFeAGga9pGjRURwoZ38OJNCHT4JLlTUB54DsyleXJXpURHhiX6Y6EEEIIIYQQ4ulA63RGZ9vdu4fY3jcBwKGRJMua44XAt8EPfAEaq8OFwHcgkWFBXYz6qjAtNVEOlgS++4ZKl6fm4iUIIYQQQgghhJihis74/mhzD43VYUwDfvJYL1Asaa4vC3wjjCXV5YcTGdrrYwCsbImXZXz7xlOFy/cPlZdACyGEEEIIIYSYnyo68N19ZIpLlzextr2OJ/2mVUHGt64k8C0vdS4Gvkubqzk0WhL4jqUBuGpNCwdGkjiuNyevQwghhBBCCCHEmavYwNd2XA6NJFnRGmdFS7xQshxMW9RSEy3cNih19jyPgUSGhXUq8G2uUZlgz1MBbt94muZ4hPMX1pGzXfrH03P8qoQQQgghhBBCnK6KDXx7x9LkHY9VLTWsaInTM5Ym77iMTGUBaIwXM75N8QhjqTyJdJ503ilkfJuqI+Qcl6TfFbp3LM2ixira6lTQPOw/1qmMTGXLxgoLIYQQQgghhJg7FRv4BmNzV7bGWd4Sx3E9ekZTjCRz1MZCRENW4bYN1RFytltoXhUEvoU5fv0scd94mkUNVTTFVeAbZI9P5a+/u4Xn/dfd/PftXbPz4oQQQgghhBBCTFvFBr5B1+WVrSrjCyoYHknmaI5Hym7bWK2yv9t6xwFY3FgNQIN/+XhKlTv3B4GvHxBPJ/BN5xweOTAKwGPdYzN9WfPeu7/3OP/1OwnwhRBCCCGEEPNHxU5n1D2aoi4WoqkkyD0wnGQ0maW5ZHwvqIwvUGiAtaSxCqBw39FUjsmsTSbv0lYXpanGzwSnTh34Pt6jgt14xKL/HJj798necSYyed2rIYQQQgghhBAFFZvxPTKRKSlZDlNfFVYZ36lcWTAMxQB3e1+C6ohVOB8ExOOpHJMZG4C6WJh4xCJimYxMI+O7+eAYhgHXX9BB33gat8I7QWdtl1TW0b0aQgghhBBCCFFQsYHv0GSW1lqV2TUMg+UtcQ6OqFLnlprjlzrvPjLFksZqDMMou3wsmWPKD3xrY2EMw6AxHi6M/QXV6OqurkEmj8p27hyYYHlznHV+J+jpBMtPZ5m8QzJn614NIYQQQgghhCio3MB3KktbbaxwfmVLnP1DSUaTOZrjxy91BljslzkD1FeFMQwYTeULAW1tTFWHN8WjjCbVZYlUnhd/+h7e+PVHuem2XWWPPTKVo602SkeDetzSKZCCaZIqSdZ2SeUk4yuEEEIIIYSYPyoy8PU8j8GJYsYXYHlznIFEBsf1jil1DppYASxpqi4shyyTuli4rNS5phD4hhlNqumMPnfXHkaSORY1VPHTLX1MZYsZz+FklpaaKIv8wLfPD3x/9ngfmz5ye+F8JfA8T2V8s5LxFUIIIYQQQswfFRn4TmZtsrZLW0ngu6I1XlhuPqrUOWyZPH9tGwCr22rKrmusDjOWyhcaNtX5gW9jtZr7F+APuwZ53nmt/O+rLyKZc/jOQ4cK9x+ZytFcEykEvv3jaTJ5h5tu28VkxubL9+yfrZetne16NHtjhHIJ3asihBBCCCGEEAUVGfgOTqhMbGnGd0VzMfA9f2HdMff56usv4bb3PptXXrKk7PLGeKQs41sbU9nh5niE0WQOz/MYSGRY3hLnoiUNvGDdAj59x256x1LkbJdEOk9LTZS6qhDxiEXfeJrf7zjCQCJD54Javv9oN1m7MkqDM3mHL0Y+zfvcb1RkGbcQZ+qRA6P88NEe3ashhBBCCHHOqsjAd2jy2MD3vPYarj1/AV97wyWsWVB7zH0Mw2Btex2RUPm/pLFaBbhB+XIwxrcxHiGRzjOazJHKOXTUV2EYBh992Xpsx+ObDxwszPPbXBPBMAw6GqroH09z/95hamMh3n3NGjJ5l50Dk2fl/zDXMnmXJiZoIUEm7+peHSHmjQ/eciu3/vT7uldDCCGEEOKcVZHz+A5OqvlyS0udoyGLW153yWk/VkN1mK7Dk0xm8limQVXYAqC9TjXOerx7HICFDep8R0MV165fwI8e6+XFG9oBaPHnDV7
"text/plain": [
"<Figure size 1152x648 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(16,9))\n",
"lfp = data_loader.lfp('1833-200619-2', 6)\n",
"# lfp = data_loader.lfp('1834-220319-3', 6)\n",
"# lfp = data_loader.lfp('1849-010319-4', 6)\n",
"times = np.arange(lfp.shape[0]) / lfp.sampling_rate.magnitude\n",
"clean_lfp, _ = compute_clean_lfp(lfp.magnitude[:, 0], threshold=2)\n",
"plt.plot(times,lfp[:,0])\n",
"plt.plot(times,clean_lfp)\n",
"\n",
"plt.figure(figsize=(16,9))\n",
"plt.psd(lfp[:,0].ravel(), Fs=1000, NFFT=10000)\n",
"plt.psd(clean_lfp, Fs=1000, NFFT=10000)\n",
"plt.xlim(0,100)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# plt.figure(figsize=(16,9))\n",
"\n",
"# plt.plot(times,lfp[:,0])\n",
"# # plt.plot(clean_lfp*100)\n",
"# plt.plot(times[:-1], np.diff(lfp[:,0].magnitude.ravel()))\n",
"# plt.xlim(64.5,65.5)\n",
"# # plt.ylim(-250,250)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"drive_0_channel_groups = [0, 1, 2, 3]\n",
"drive_1_channel_groups = [4, 5, 6, 7]"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4\n",
"5\n",
"6\n",
"7\n",
"4\n",
"5\n",
"6\n",
"7\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f43cc322940>]"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD3CAYAAAD7VehMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xT9f7H8ddJV7pbyigFCrK+IHuooKDgVtzee73Xva6o9zrv7w7HvY7rVe91X73uPfFeLzi4OFBBZYns/YUyCi0ddM+kaZLfHyeBUFpImzRJ28/z8fBhm3OSfBtOzvt85zHcbjdCCCG6Jku4CyCEECJ8JASEEKILkxAQQoguTEJACCG6MAkBIYTowqLDXQCv4447zt2nT59wF0MIITqMjRs3lmitewTyGhETAn369GH27NnhLoYQQnQYSqncQF9DmoOEEKILkxAQQoguTEJACCG6MAkBIYTowiQEhBCiC5MQEEKILkxCQAghujAJASGE6MIkBIQQogs74oxhpZQFeB4YA9iB67XWOT7b7wB+6fl1ntb6AaWUAeQB2zyPL9Va3xXUkgshhAiYP8tGXABYtdaTlVKTgCeA8wGUUgOBy4DjABewSCk1B6gDVmmtz22fYgshhAgGf5qDpgBfAGitlwETfbbtAc7UWju11m4gBrABE4A+SqkFSql5SikV5HILIYQIAn9CIAWo9PndqZSKBtBaO7TWJUopQyn1OLBaa70VKAAe0VpPBx4G3g12wYUQQgTOnxCoApJ9n6O1bvT+opSyAu959rnZ8/AK4BMArfUiIMvTTyCEECKC+BMCi4GzATx9Auu9Gzwn9k+AtVrrmVprp2fTfcDtnn3GAHs8zUVCCCEiiD8dw3OA05RSSwADuEYpdSeQA0QBJwFxSqmzPPvfBTwKvKuUmgE0AlcHu+BCCCECd8QQ0Fq7gBubPLzF52drC0+d0dZCCSGECA2ZLCaEEF2YhIAQQnRhEgJCCNGFSQgIIUQXJiEghBBdmISAEEJ0YRICQgjRhUkICCFEFyYhIIQQXZiEgBBCdGESAkII0YVJCAghRBcmISCEEF2YhIAQQnRhEgJCCNGFSQgIIUQXJiEghBBdmISAiFhVNgf/WpCDzeE88s5CiDaREBARa8GWYh77UvPqDzvCXRQhOi0JARGx8srrAXjxux2U1NjDXBohOicJARGx8srrscZYqHc4+ec328JdHCE6JQkBEbHyK+oZ0jOZXx3bj/d/3M2OfTXhLpIQnY6EgIhYeyvq6ZMWz22nDCUu2sI/vtDhLpIQnY6EgIhIbreb/PJ6stLi6ZEcx8yTBvHFxkJW7CoLd9GE6FQkBEREKq9zUO9w0ic9HoDrpx5Fz+Q4Hp63GbfbHebSCdF5SAiIiJTvGRnUJ80MgYTYaG4/dSirdlewMrc8nEUTolOREBARKb+iDoC+npoAwPj+aQAUV8twUSGCRUJARKS8JjUBgKS4aABqbI1hKZMQnZGEgIhI+RX1JMRGkZYQs/+xZKv5c5XNEa5iCdHpSAiIiOQdHmoYxv7H9tcE7FITECJYJARERMqvqN8/MsgrymKQEBslzUFCBFH0kXZQSlmA54ExgB24Xmud47P9DuCXnl/naa0fUErFA+8CPYFq4Cqt9b5gF150Xvnl9Yzum3bI48nWaKolBIQIGn9qAhcAVq31ZOBPwBPeDUqpgcBlwPHAJOB0pdRo4CZgvdZ6KvA2cG+wCy46r7qGRsrrHAd1CnslxUVLc5AQQeRPCEwBvgDQWi8DJvps2wOcqbV2aq3dQAxg830O8DlwatBKLDo97xyBvunNhIA1hmoJASGCxp8QSAEqfX53KqWiAbTWDq11iVLKUEo9DqzWWm9t8pxqIDWYhRadW17FocNDvZLjoqmR0UFCBI0/IVAFJPs+R2u9/1JMKWUF3vPsc3Mzz0kGKgIvqugq9s8WbqYmIH0CQgSXPyGwGDgbQCk1CVjv3aCUMoBPgLVa65laa2fT5wBnAT8ErcSi08uvqCfaYtAz2XrINukTECK4jjg6CJgDnKaUWgIYwDVKqTuBHCAKOAmIU0qd5dn/LuAF4C2l1CKgAbg06CUXndbeinp6p1mJshiHbEuyRssQUSGC6IghoLV2ATc2eXiLz8+HXq6Zft7WQomuLb+8vtn+ADBnDdc0NOJyubE0ExJCiNaRyWIi4uRXmPcRaE5yXDRuN9Q2SG1AiGCQEBARxeF0UVRlo28LIZBklaUjhAgmCQERUQorbbjczY8MggPrB8kIISGCQ0JARJQDS0gnNLs92SohIEQwSQiIiJJf0fIcATgQAtIcJERwSAiIiLLXEwK9U5sfdJYUZ95ToFpmDQsRFBICIqLkl9fTIzkOa0xUs9v31wSkOUiIoJAQEBElv6LlOQIgo4OECDYJARFRmruZjK/EWOkYFiKYJARExHC53EesCURZDJLiZBE5IYJFQkBEjJJaOw2NrsOGAHgXkZOOYSGCQUJARIz9S0gfKQSsspKoEMEiISAixt4KG9DyHAEvaQ4SIngkBETEKK42Q6BXSksL05rkxjJCBI+EgIgY5bUNWAxIjY857H7J0hwkRNBICIiIUVbXQFpCbLM3k/FlNgdJx7AQwSAhICJGeZ2DtITD1wLAc2MZaQ4SIigkBETEKK9toFtC7BH3S4qLprbBidPlDkGphOjcJARExCirNZuDjkRWEhUieCQERMSoqHPQLdGf5iAJASGCRUJARAS3201ZXQPpif40B5lBIf0CQgROQkBEhLoGJw2NLtL96RPYf3cxGSEkRKAkBEREKK9rAPC7YxigWpqDhAiYhICICOW15lW9P81BKXJjGSGCRkJARIQyT00g3Y95AnJjGSGCR0JARIQKbwj41TEsfQJCBIuEgIgIZbX+9wkkxkZjGNIcJEQwSAiIiFBe24BhQMoRFo8DsFgMkmKjpWNYiCCQEBARobzOQVp8zBEXj/NKkuWkhQgKCQEREfydKOaVbI2W5iAhgkBCQESEiroGvyaKeZn3GZYQECJQ0UfaQSllAZ4HxgB24HqtdU6TfXoAi4HRWmubUsoA8oBtnl2Waq3vCmrJRadSVus44r2FfSVZY6j0jCgSQrTdEUMAuACwaq0nK6UmAU8A53s3KqXOAB4FMn2eMwhYpbU+N5iFFZ1XeW0Do/qk+L1/clw0eeV17VgiIboGf5qDpgBfAGitlwETm2x3AacCZT6PTQD6KKUWKKXmKaVUMAorOie32015K5uDpE9AiODwJwRSgEqf351Kqf01CK31fK11aZPnFACPaK2nAw8D7wZcUtFp1Tuc2BtdreoYNm8xKSEgRKD8CYEqINn3OVrrI337VgCfAGitFwFZnn4CIQ7RmoliXknWaOodThqdrvYqlhBdgj8hsBg4G8DTJ7Dej+fcB9zuec4YYI/WWu4FKJpVUWcu/+DP/YW9kq3mvrV2Z7uUSYiuwp+O4TnAaUqpJYABXKOUuhPI0Vp/2sJzHgXeVUrNABqBq4NRWNE57a8JtGaewP7lpB2ktiI8hBAHO2IIaK1dwI1NHt7SzH4DfH4uB2YEWjjRNXjvJeDP/YW9DtxYRvoFhAiETBYTYVfelpqALCctRFBICIiwK6tzYBiQ6sficV7e5aRlmKgQgZEQEGFXXttAaisWj4MDNYEquaeAEAGREBBhV17X0KrhoQBJcWatQZqDhAiMhIAIu/JWriAKPn0C0hwkREAkBETYldU6/Lq3sK+E2CgMQ0YHCREoCQERdq1dRhrAMAxZTlqIIJAQEGFXVtvQquGhXinWGKkJCBEgCQERVvUN5uJxrZko5mUuIiejg4QIhISACKuyOu9EsdYv/ZBkleYgIQIlISDCyjtbuLV9AuC5p4CEgBABkRAQYeVdN6i1Q0RB7ikgRDBICIiwKguwJiAhIERgJAREWLVl8Tgvc4iodAwLEQgJARFW5W1YPM4r2RqDzeHCIXcXE6LNJAREWJXXtX7xOC9ZSVSIwEkIiLAqq2394nFeSXJPASECJiEgwqqiztGqewv7SpG7iwkRMAkBEVZtXTICDiwnLbO
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAEMCAYAAAAI6znIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydeXwTdf7/X0maNE2apkd60oNeM21BsBxSKOCBAh7ghbC7ghd4oSvuIq7Xuofu6uKBKLriKsqiP8Fr+XosKyq6oMDKImfTTu/0PtM2adqc8/79MSS20JYkkx64eT4en0ePzHzmPZOZ93w+7+sjISIECRIkiBikoy1AkCBBzn2CiiRIkCCiCSqSIEGCiCaoSIIECSKaoCIJEiSIaIKKJEiQIKIJGW0Bgow9WJYNBfAmgAwAJgD3AIgBsBGAE8BujuP+wLKsFsC2U7ut4DiuazTkDTL6BBVJkIG4HUA3x3EFLMuyADYBiAdwPYBKAJ+xLJsPIBPAMxBGtpcC+HCU5A0yygSnNkEGIg/ALgDgOI4DMB1AKMdxFRzHEYDPISiOLwA8AODXAL4cJVmDjAGCI5IgA3EUwFUsy+4EMAOAFkBFn8/NADJOTWWuHgX5gowxgiOSIAOxBYJtZB+AawEcA6Du87kGQOcoyBVkjBJUJEEGYjqArziOmw3gfQClAOwsy2ayLCsBsACCkgkSBAAgCSbtBTkdlmV1ALZDGIV0AlgJIBXACwBkELw2j46ehEHGGkFFEiRIENEEpzZBggQRTVCRBAkSRDRBRRIkSBDRBBVJkCBBRBNUJEGCBBFNUJEECRJENEFFEiRIENEEFUmQIEFEE1QkQYIEEU0w+zcIAIBlWQWABACJABJDQkKSVSpVjFQqVUilUoVMJlNIJBI5ABkROYjIwfO8zeVy2Z1Op9VisdQDaOzT2k+VHAjyP0AwRP5/hFPVzKYoFIpJWq02OyQkJJ3n+SSe52NkMlloaGioPDY2FnFxcSE6nS5Up9PJw8LCJBKJBFKp1NMcDgdCQkLA83y/1tXVxbe1tdlaWlrsra2trs7OTjidThvP8z0ymawZQK3Vaq3s6uoqJqIjAEo4jnOO8mUJEiCCiuQnCMuykQCmhIeHzwoPD7/Y6XTmREZGKidOnKjIyMhQa7VaSUREBKKiohAVFYXQ0FDI5XKEhIRALpdDJpNBJpNBKpVCJpNBIpF4WllZGbKysuC+b3ieh8vlgsvlAs/zcDqdcDgcnp8WiwVGoxGdnZ0wm81obW116fV6c0VFhcvhcLRLpdJDra2tX/E8/18AxUHlcm4SVCTnOCzLygAUhIeHzzulNNjo6OjQvLw8ZWZmpiotLQ2pqalQq9VQKpWeJpfLIZFIfD5ecXExcnNzfd6P53nYbDZYrVZPMxqNqKysRE1NjaukpMRSXl7ucDgcHVKp9FB7e/uXTqdzN8dxdT4fLMiIE1Qk5yAsy2okEsn8uLi4W4nogilTpignT56sSUlJwbhx4xAREQG1Wg2VSgWFQuGXwhgMfxXJYPA8j97eXvT09KC7uxtGoxE1NTUwGAyuffv2dTU3N3c6nc5/dHR0vAvgh6DdZWwSVCTnCCzLpqhUquvCw8NvUiqV42fPnh2en5+vyM7ORmRkJDQaDVQqFaTS4XXEBVqRDITD4UB3dzdMJhNaWlpw8uRJHDx40Hzy5MlemUz27+bm5rcA7OE4zjqsggTxmqAiGcOwLHt+VFTUcplMdl1iYmLknDlztJMnT5amp6dDq9VCo9FAJpONqEwjoUhOx2azoaurC+3t7dDr9Th8+HDvgQMHelwuV0lHR8dWu93+IcdxxhEVKkg/gopkjMGyrFalUt2sUqnWMAwTffHFF0fm5eUhISEBUVFRUKlUoqcqdrsdzc3NaGxsRFNTExobG9HY2IiGhoZ+v3d3d+OCCy7ApZdeinnz5mHKlCkoLS0dcUXSF5fLha6uLnR0dKCqqgrHjh1z7tq1q8tmsx1qaWl5GsDe4PRn5AkqkjHAqTqoF8TFxT0kl8vnXHnlldqZM2eGJCcnIzo6Gmq1elDlQUQoLS3F7t278eWXX6K1tdXjqu3ruuV5Hq2trWhqakJ7e/sZ/UgkEsTGxiIxMRFJSUlITExEaGgovv32W5w4cQIAEBkZialTp+Laa6/FkiVLEB8fP5yX5ay4XC50dnaira0NHMdh9+7d5iNHjnTa7fbXTCbTqxzHtY2qgP9DBBXJKMKyrCI0NPRnGo3m0dzc3IQrrrgiYuLEidDpdIiMjBzU3mE0GvHVV19h9+7d2L17N2pqagAAmZmZSE9PBxF54jvcvwOATqdDYmIiEhISzvgZHx8PuVw+4PGam5uxZ88efPXVV/jXv/6F+vp6REZG4sUXX8Ty5csDasz1F5vNhvb2dtTX1+M///mPa+fOnZ1Wq/Wb1tbWP3Acd2K05fupE1QkowDLsrGRkZG/DgkJuXX+/PmRl1xySej48eMRGxsLpVI55L4fffQRfv7zn8Nut0Or1WLevHmYP38+LrvsMmRkZAy77MXFxXC5XLjzzjuxf/9+LF68GJs3b0ZCQsKwH9sbiAhdXV1oaWnBsWPHsHPnzk6DwWAwGo2/d7lcH3Mcx4+2jD9JiCjYRqgxDBMxY8aM5y+++OK2zZs3O3/44Qdqbm4mp9NJ3rBz504KCQmhmTNn0v79+8nhcHi1XyDR6/VEROR0Oum5554jpVJJ0dHR9M477xDP8yMuz1D09PRQdXU1ffXVV7Ru3brumTNnVuTl5V3JMIyExsD98FNqoy7A/0JjGCZ02rRpD8+ePbtl06ZNzhMnTlB7e7tPD94nn3xCcrmcZsyYQV1dXV7vF2jcisRNSUkJFRQUEAC65pprqKmpaZQkGxy73U61tbX0zTff0OrVq82zZs06xjBMAY2Be+On0kZdgJ9yYxhGNnny5JWzZs1qeOqpp+xHjx6ljo4On9/cu3btIoVCQdOmTaOOjg6f9g00pysSImF0sn79egoNDaWUlJQBtxkLOBwOqquro127dtHy5ctNs2bN+jfDMLk0Bu6Vc70FbSTDAMuyErlcvlir1b5w0UUXJV599dWh7sAxXw2TX3zxBRYtWoS8vDx8+eWXiI6OHiapvWOoOJIjR47g8ssvh91ux6effopZs2aNsHTe4XQ60dTUhO+//x5///vfu1pbW/e0tbWt4TiudrRlO1cJKpIAw7JsoU6n2zx16tTxS5YsUefm5kKn0/nl2SgqKsK0adPAMAz27NmDmJiYYZDYN84WkFZZWYkFCxagvr4eO3bswKJFi0ZQOt+w2+2oq6vDwYMH8dZbb3V0d3e/397e/gjHcWf6x4MMSVCRBAiWZbUxMTGvp6enX3bzzTdrJ0yYgISEBFGRpzfeeCM+/vhjlJeXj1jMBhH1y+h1Z/W6P6upqUFqaqonG9idIdw3Y7ilpQVXXnkljhw5gs2bN2PlypUjIru/9PT0oLq6GgcOHOC3bt3aZjKZ7j969Oi7oy3XuURQkQSACRMmLIiKinrz1ltvjZ8zZ440NTX1rG7cs1FVVYXs7Gzcf//9ePbZZwMkqaAMHA6HJwPXbrfDZrPB4XB4FIZMJkNISIhHSbiD2yQSCVpbW6HT6Txz48GUjs1mw/3334+9e/fi4YcfxiOPPIKwsLARD+n3FiKC0WhERUUF3nrrLcvRo0f3t7e338hxXOtoy3YuEFQkImBZNiImJub1rKyshXfccYdmwoQJftlBBuLee+/Fa6+9hqqqKowbN86vPogIVqsVFosFPT096O3thcvlgkKhgFKpRGhoKEJDQ6FQKKBQKLxK+PMm14aI4HK5YLFYcOedd2LHjh248847sWbNGhARFAoFVCoV1Go11Go1QkLGTqE+l8uF+vp67N27lzZv3txmNpvvO3r06PbRlmusc04qEpZlZwD4C8dxF7EsmwfgNQASAGUAVnEc52RZdiOA2QDMp3a7GkAGgA0AjnIcd78YGSZMmDA/Kipq68033xx/2WWXSZKTkwP2QLS2tiItLQ0/+9nPsGXLFq/3IyL09vbCZDLBbDbD4XB
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# action_id_0, channel_0, unit_0 = '1833-200619-1', 6, 163\n",
"# action_id_1, channel_1, unit_1 = '1833-200619-2', 6, 28\n",
"action_id_0, channel_0, unit_0 = '1834-220319-3', 2, 46\n",
"action_id_1, channel_1, unit_1 = '1834-220319-4', 2, 60\n",
"\n",
"# change data loader to get all LFPs and then selecte the best form the other\n",
"lfp_0 = data_loader.lfp(action_id_0, channel_0)\n",
"lfp_1 = data_loader.lfp(action_id_1, channel_1)\n",
"\n",
"# select best channel among other drive\n",
"if channel_0 in drive_0_channel_groups:\n",
" lfps = []\n",
" for ch in drive_1_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id_0, ch))\n",
"elif channel_0 in drive_1_channel_groups:\n",
" lfps = []\n",
" for ch in drive_0_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id_0, ch))\n",
" \n",
"# merge lfp of other drive into a single AnalogSignal\n",
"lfp_arrays = np.hstack(lfps).as_array()\n",
"lfp_0 = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,\n",
" t_start=lfps[0].t_start)\n",
"\n",
"# select best channel among other drive\n",
"if channel_1 in drive_0_channel_groups:\n",
" lfps = []\n",
" for ch in drive_1_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id_1, ch))\n",
"elif channel_1 in drive_1_channel_groups:\n",
" lfps = []\n",
" for ch in drive_0_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id_1, ch))\n",
"# merge lfp of other tetrodes into a signle AnalogSignal\n",
"lfp_arrays = np.hstack(lfps).as_array()\n",
"lfp_1 = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,\n",
" t_start=lfps[0].t_start)\n",
"\n",
"\n",
"sample_rate_0 = lfp_0.sampling_rate\n",
"sample_rate_1 = lfp_1.sampling_rate\n",
"\n",
"lim_0 = get_lim(action_id_0)\n",
"lim_1 = get_lim(action_id_1)\n",
"\n",
"sptrs_0 = data_loader.spike_trains(action_id_0, channel_0)\n",
"\n",
"sptrs_1 = data_loader.spike_trains(action_id_1, channel_1)\n",
"\n",
"cleaned_lfps_0, sptr_0, best_channel_0 = prepare_spike_lfp(lfp_0, sptrs_0[unit_0], *lim_0)\n",
"\n",
"cleaned_lfps_1, sptr_1, best_channel_1 = prepare_spike_lfp(lfp_1, sptrs_1[unit_1], *lim_1)\n",
"\n",
"coher_0, freq_0 = compute_spike_lfp_coherence(cleaned_lfps_0, sptr_0, 4096)\n",
"\n",
"coher_1, freq_1 = compute_spike_lfp_coherence(cleaned_lfps_1, sptr_1, 4096)\n",
"\n",
"spike_phase_0, filtered_lfp_0 = compute_spike_phase(cleaned_lfps_0, sptrs_0[unit_0], flim=[6,10])\n",
"\n",
"spike_phase_1, filtered_lfp_1 = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[6,10])\n",
"\n",
"# spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[10.5,11.5])\n",
"spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[29.5,30.5])\n",
"\n",
"plt.figure()\n",
"plt.plot(freq_0, coher_0.ravel())\n",
"plt.plot(freq_1, coher_1.ravel())\n",
"plt.xlim(0,20)\n",
"\n",
"plt.figure()\n",
"bins_0, kde_0 = vonmises_kde(spike_phase_0, 100)\n",
"ang_0, vec_len_0 = spike_phase_score(bins_0, kde_0)\n",
"plt.polar(bins_0, kde_0, color='b')\n",
"plt.polar([ang_0, ang_0], [0, vec_len_0], color='b')\n",
"\n",
"bins_1, kde_1 = vonmises_kde(spike_phase_1, 100)\n",
"ang_1, vec_len_1 = spike_phase_score(bins_1, kde_1)\n",
"plt.polar(bins_1, kde_1, color='r')\n",
"plt.polar([ang_1, ang_1], [0, vec_len_1], color='r')\n",
"\n",
"bins_1_stim, kde_1_stim = vonmises_kde(spike_phase_1_stim, 100)\n",
"ang_1_stim, vec_len_1_stim = spike_phase_score(bins_1_stim, kde_1_stim)\n",
"plt.polar(bins_1_stim, kde_1_stim, color='k')\n",
"plt.polar([ang_1_stim, ang_1_stim], [0, vec_len_1_stim], color='k')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# TODO fix artefact stuff from phase precession"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"NFFT = 8192\n",
"theta_band_f1, theta_band_f2 = 6, 10 \n",
"drive_0_channel_groups = [0, 1, 2, 3]\n",
"drive_1_channel_groups = [4, 5, 6, 7]"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"coherence_data, freqency_data = {}, {}\n",
"theta_kde_data, theta_bins_data = {}, {}\n",
"stim_kde_data, stim_bins_data = {}, {}\n",
"\n",
"def process(row):\n",
" action_id = row['action']\n",
" channel_group = row['channel_group']\n",
" unit_name = row['unit_name']\n",
" \n",
" name = f'{action_id}_{channel_group}_{unit_name}'\n",
" \n",
" # select best channel among other drive\n",
" if channel_group in drive_0_channel_groups:\n",
" lfps = []\n",
" for ch in drive_1_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id, ch))\n",
" elif channel_group in drive_1_channel_groups:\n",
" lfps = []\n",
" for ch in drive_0_channel_groups:\n",
" lfps.append(data_loader.lfp(action_id, ch))\n",
" \n",
" # merge lfp of other tetrodes into a signle AnalogSignal\n",
" lfp_arrays = np.hstack(lfps).as_array()\n",
" lfp = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,\n",
" t_start=lfps[0].t_start) \n",
" sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
" \n",
" lim = get_lim(action_id)\n",
" \n",
" cleaned_lfp, sptr, best_channel = prepare_spike_lfp(lfp, sptr, *lim)\n",
" \n",
" p_xys, freq = compute_spike_lfp_coherence(cleaned_lfp, sptr, NFFT=NFFT)\n",
" \n",
" p_xy = p_xys.magnitude.ravel()\n",
" freq = freq.magnitude\n",
" \n",
" theta_f, theta_p_max = find_theta_peak(p_xy, freq, theta_band_f1, theta_band_f2)\n",
" \n",
" theta_energy = compute_energy(p_xy, freq, theta_band_f1, theta_band_f2) # theta band 6 - 10 Hz\n",
" \n",
" theta_half_f1, theta_half_f2 = compute_half_width(p_xy, freq, theta_p_max, theta_f)\n",
" \n",
" theta_half_width = theta_half_f2 - theta_half_f1\n",
" \n",
" theta_half_energy = compute_energy(p_xy, freq, theta_half_f1, theta_half_f2) # theta band 6 - 10 Hz\n",
" \n",
" theta_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[theta_band_f1, theta_band_f2])\n",
" theta_bins, theta_kde = vonmises_kde(theta_spike_phase)\n",
" theta_ang, theta_vec_len = spike_phase_score(theta_bins, theta_kde)\n",
" theta_kde_data.update({name: theta_kde})\n",
" theta_bins_data.update({name: theta_bins})\n",
"\n",
" # stim\n",
" \n",
" stim_freq = compute_stim_freq(action_id)\n",
" \n",
" stim_p_max = compute_stim_peak(p_xy, freq, stim_freq)\n",
" \n",
" stim_half_f1, stim_half_f2 = compute_half_width(p_xy, freq, stim_p_max, stim_freq)\n",
" stim_half_width = stim_half_f2 - stim_half_f1\n",
" \n",
" stim_energy = compute_energy(p_xy, freq, stim_half_f1, stim_half_f2)\n",
" \n",
" if np.isnan(stim_freq):\n",
" stim_spike_phase, stim_bins, stim_kde, stim_ang, stim_vec_len = [np.nan] * 5\n",
" else:\n",
" stim_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[stim_freq - .5, stim_freq + .5])\n",
" stim_bins, stim_kde = vonmises_kde(stim_spike_phase)\n",
" stim_ang, stim_vec_len = spike_phase_score(stim_bins, stim_kde)\n",
" stim_kde_data.update({name: stim_kde})\n",
" stim_bins_data.update({name: stim_bins})\n",
" \n",
" coherence_data.update({name: p_xy})\n",
" freqency_data.update({name: freq})\n",
" \n",
" result = pd.Series({\n",
" 'best_channel': best_channel,\n",
" 'theta_freq': theta_f,\n",
" 'theta_peak': theta_p_max,\n",
" 'theta_energy': theta_energy,\n",
" 'theta_half_f1': theta_half_f1, \n",
" 'theta_half_f2': theta_half_f2,\n",
" 'theta_half_width': theta_half_width,\n",
" 'theta_half_energy': theta_half_energy,\n",
" 'theta_ang': theta_ang, \n",
" 'theta_vec_len': theta_vec_len,\n",
" 'stim_freq': stim_freq,\n",
" 'stim_p_max': stim_p_max,\n",
" 'stim_half_f1': stim_half_f1, \n",
" 'stim_half_f2': stim_half_f2,\n",
" 'stim_half_width': stim_half_width,\n",
" 'stim_energy': stim_energy,\n",
" 'stim_ang': stim_ang, \n",
" 'stim_vec_len': stim_vec_len\n",
" })\n",
" return result"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6e6aa8f543ce4316bb63d5c3b6761864",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(IntProgress(value=0, max=1284), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1577: RuntimeWarning: invalid value encountered in true_divide\n",
" Cxy = np.abs(Pxy)**2 / Pxx / Pyy\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in true_divide\n"
]
}
],
"source": [
"results = units.merge(\n",
" units.progress_apply(process, axis=1), \n",
" left_index=True, right_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"# coher, freqs = {}, {}\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",
" \n",
"# name = f'{action_id}_{channel_group}_{unit_name}'\n",
" \n",
"# lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response\n",
" \n",
"# sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
" \n",
"# lim = get_lim(action_id)\n",
"\n",
"# p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)\n",
" \n",
"# snls = signaltonoise(clean_lfp)\n",
"# best_channel = np.argmax(snls)\n",
"# snl = snls[best_channel]\n",
"# p_xy = p_xys[:,best_channel].magnitude\n",
"# freq = freq.magnitude\n",
" \n",
"# coher.update({name: p_xy})\n",
"# freqs.update({name: freq})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pd.DataFrame(coherence_data).to_feather(output / 'data' / 'coherence.feather')\n",
"pd.DataFrame(freqency_data).to_feather(output / 'data' / 'freqs.feather')\n",
"pd.DataFrame(theta_kde_data).to_feather(output / 'data' / 'theta_kde.feather')\n",
"pd.DataFrame(theta_bins_data).to_feather(output / 'data' / 'theta_bins.feather')\n",
"pd.DataFrame(stim_kde_data).to_feather(output / 'data' / 'stim_kde.feather')\n",
"pd.DataFrame(stim_bins_data).to_feather(output / 'data' / 'stim_bins.feather')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Save to expipe"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"action = project.require_action(\"stimulus-spike-lfp-response-other-drive\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"action.modules['parameters'] = {\n",
" 'NFFT': NFFT,\n",
" 'theta_band_f1': theta_band_f1,\n",
" 'theta_band_f2': theta_band_f2\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"action.data['results'] = 'results.csv'\n",
"results.to_csv(action.data_path('results'), index=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"copy_tree(output, str(action.data_path()))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"septum_mec.analysis.registration.store_notebook(action, \"10-calculate-stimulus-spike-lfp-response-other-drive.ipynb\")"
]
}
],
"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": 4
}