Mikkel
/
septum-mec
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
septum-mec/actions/stimulus-lfp-response-mec-n.../data/20_stimulus-lfp-response-me...

1416 lines
259 KiB
Plaintext
Raw Permalink Normal View History

actions/stimulus-lfp-response-mec-no-zscore/
2021-03-10 12:32:52 +00:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import expipe\n",
"import pathlib\n",
"import numpy as np\n",
"import spatial_maps.stats as stats\n",
"import septum_mec\n",
"import septum_mec.analysis.data_processing as dp\n",
"import septum_mec.analysis.registration\n",
"import head_direction.head as head\n",
"import spatial_maps as sp\n",
"import speed_cells.speed as spd\n",
"import re\n",
"import joblib\n",
"import multiprocessing\n",
"import shutil\n",
"import psutil\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib\n",
"from distutils.dir_util import copy_tree\n",
"from neo import SpikeTrain\n",
"import scipy\n",
"import seaborn as sns\n",
"from tqdm.notebook import tqdm_notebook as tqdm\n",
"tqdm.pandas()\n",
"\n",
"from spike_statistics.core import permutation_resampling_test\n",
"\n",
"from spikewaveform.core import calculate_waveform_features_from_template, cluster_waveform_features\n",
"\n",
"from septum_mec.analysis.plotting import violinplot, despine"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#############################\n",
"\n",
"perform_zscore = False\n",
"\n",
"if not perform_zscore:\n",
" zscore_str = \"-no-zscore\"\n",
"else:\n",
" zscore_str = \"\"\n",
"\n",
"#################################"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"plt.rc('axes', titlesize=12)\n",
"plt.rcParams.update({\n",
" 'font.size': 12, \n",
" 'figure.figsize': (6, 4), \n",
" 'figure.dpi': 150\n",
"})\n",
"\n",
"output_path = pathlib.Path(\"output\") / (\"stimulus-lfp-response-mec\" + zscore_str)\n",
"(output_path / \"statistics\").mkdir(exist_ok=True, parents=True)\n",
"(output_path / \"figures\").mkdir(exist_ok=True, parents=True)\n",
"output_path.mkdir(exist_ok=True)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"data_loader = dp.Data()\n",
"actions = data_loader.actions\n",
"project = data_loader.project"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"identify_neurons = actions['identify-neurons']\n",
"sessions = pd.read_csv(identify_neurons.data_path('sessions'))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"lfp_action = actions['stimulus-lfp-response' + zscore_str]\n",
"lfp_results = pd.read_csv(lfp_action.data_path('results'))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"lfp_results = pd.merge(sessions, lfp_results, how='left')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"lfp_results = lfp_results.query('stim_location!=\"ms\"')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def action_group(row):\n",
" a = int(row.channel_group in [0,1,2,3])\n",
" return f'{row.action}-{a}'\n",
"lfp_results['action_side_a'] = lfp_results.apply(action_group, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"lfp_results['stim_strength'] = lfp_results['stim_p_max'] / lfp_results['theta_energy']"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>action_side_a</th>\n",
" <th>channel_group</th>\n",
" <th>signal_to_noise</th>\n",
" <th>stim_strength</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>71</th>\n",
" <td>1833-010719-1-0</td>\n",
" <td>7</td>\n",
" <td>0.001902</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>67</th>\n",
" <td>1833-010719-1-1</td>\n",
" <td>3</td>\n",
" <td>0.003522</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>583</th>\n",
" <td>1833-020719-1-0</td>\n",
" <td>7</td>\n",
" <td>-0.002942</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>579</th>\n",
" <td>1833-020719-1-1</td>\n",
" <td>3</td>\n",
" <td>0.012323</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>375</th>\n",
" <td>1833-020719-3-0</td>\n",
" <td>7</td>\n",
" <td>-0.002042</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" action_side_a channel_group signal_to_noise stim_strength\n",
"71 1833-010719-1-0 7 0.001902 NaN\n",
"67 1833-010719-1-1 3 0.003522 NaN\n",
"583 1833-020719-1-0 7 -0.002942 NaN\n",
"579 1833-020719-1-1 3 0.012323 NaN\n",
"375 1833-020719-3-0 7 -0.002042 NaN"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# lfp_results_hemisphere = lfp_results.sort_values(\n",
"# by=['action_side_a', 'stim_strength', 'signal_to_noise'], ascending=[True, False, False]\n",
"lfp_results_hemisphere = lfp_results.sort_values(\n",
" by=['action_side_a', 'channel_group'], ascending=[True, False]\n",
").drop_duplicates(subset='action_side_a', keep='first')\n",
"lfp_results_hemisphere.loc[:,['action_side_a','channel_group', 'signal_to_noise', 'stim_strength']].head()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"colors = ['#1b9e77','#d95f02','#7570b3','#e7298a']\n",
"labels = ['Baseline I', '11 Hz', 'Baseline II', '30 Hz']\n",
"# Hz11 means that the baseline session was indeed before an 11 Hz session\n",
"queries = ['baseline and i and Hz11', 'frequency==11', 'baseline and ii and Hz30', 'frequency==30']"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# prepare pairwise comparison: same animal same side same date different sessions"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"def make_entity_date_side(row):\n",
" s = row.action_side_a.split('-')\n",
" del s[2]\n",
" return '-'.join(s)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"lfp_results_hemisphere['entity_date_side'] = lfp_results_hemisphere.apply(make_entity_date_side, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/lib/histograms.py:898: RuntimeWarning: invalid value encountered in true_divide\n",
" return n/db/n.sum(), bin_edges\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFICAYAAACSp82YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5xcVf3/8dembbIJIY0mAQIEPoQuHYFfgl+a0sQIX5AAAQFBAUGkKFWk61eUIigKQUAQkBqKIk0BQaSH8qEmIJ00SDbZJLv7++OcYW9m78zO7s7uzN19Px+PfdyZe88998xMMp85555S09zcjIiIiFS3PpUugIiIiLRNAVtERCQDFLBFREQyQAFbREQkAxSwRUREMkABW0REJAMUsEVERDJAAVtERCQDFLBFREQyQAFbREQkAxSwRUREMkABW0REJAMUsEVERDJAAVtERCQDFLClapnZhin7pphZc/xbsRLlyitPqzKKiHSFfpUugEg+M1sWOAv4PlX6bzT+WPgFsA2weoWLIyK9QFV+GUqv90vgkEoXog3XA18FZlS6ICLSO6hJXKpR30oXoARZKKOI9CAK2CIiIhmggC0iIpIBNc3NzZUugwgAZnYmcEaBw4+4+wQzmwIcFPetBAwETgB2AVYGPgNeBqYA17h7wX/gZtYfmAzsDWwIjADmAM8DtwBXu/uivHOS1893jbtPzku/AnAY8D+AxWssAWYCTwE3An8pVs5yMLNxwFGxHKOBGuBd4CHgEnd/ucB5Uwiv93l339jM1gKOBXYmvN/1wAvAHwmvv6lIGdr9fsfzxgBvx6d7AR8Q+jlsAiwgfN7fc/cXEudsCBxH6GewUrzOk8DF7n6/md0XX8Mj7j4hnnMpoaMjwFfd/aEir2Uf4M/x6e7uPrVQWpFyUQ1bsmwi8ArwPWANoBZYDhgPXA3cbWapHSvNbCwhUPwO2BFYAegfz98BuAJ43syso4Uzs4MIgeZnwARC4KgFBgOrxvLfDNxhZl12T9zMTgNeJLxPFq9fFx8fAbxoZmeaWU0b+UwEnov5rEn4sTSC8NquAv5qZrUFzi3X+/1lwo+Mr8TrDwc2At5KXOtg4BnCj4NVE9fZDfibmf28QN5TEo/3b6McB8TtR8B9baQVKQsFbKkmVxC+kO9K7Pty/Ds0Jf2lhJriZYQa9vbAKcD8ePxrwA/zT4pDsv4JjAMWxfN3BbYA9gSuARqBdYCHzGylxOmnx/I8HZ9/kCjj6YlrfJUQAAYBswjD1L4GbEUI1JfFawPsXuD1dVpstTiL0EnuBUKA/gqwLfAD4E3C98AZFG7dgBD4ro+PLyTU1LcFTia0akAIvMemlKEz73e+Uwmf+Y8JQ+oOBM5093nxWhMJPx76xnKdAWxHqGn/X7z+jwj/Vpbi7v8BpsWn3yry42M5wr83gOvdfUmR8oqUjYZ1SdVw9w+BD81sVmLfc0VOaSA0XT6e2PewmT0I/Cs+P4gQYJKuAFYkfKHv4O5P5R2/08xuAe4k1IovAvaN5XkHeMfM5sW0iwqU8ay4XQLsHINB0q2xWTb342Rv4LdFXmu7mdkmwGnx6bXAIXnB5TEz+wMwlVBLPt3MbirQPD6c8ENo27zX+5iZPQI8Tgikk4EL8s7t8Pudog9wjLtfGZ9/8dmb2WDgV/HprFjWVxLnPmRmtwP3E2rnaaYQxtcvS6iR/yUlzb60fHdeUyAfkbJTDVuy7NK8YA2Auz9BaBIFWCfZLG5mawN7xKfnpASPXB5Tafky3tvMvlRqocysjtD0PQu4KyVYJ68xJz5dudT82+F4wv/xmcARaTVBd59PGPPeTAi4RxfJ7zdpP07i+527f7xOsmbaBe/3AgoHyYmE+/MAP84L1rnrPAqcU+B8gOsIP7KgcLN4rjn8ueR9c5GupoAtWXZPkWOvx20fQm0p5+uEwAShplVK/n0INdCSuHu9u2/u7iOBb7WR/MO4TW1+7ah4P/pr8elj7l5fKK27v03oCwChqbuQvxU59mbi8ZDE43K/38+kdUyLdo/bxbQ036f5XaED7v4RcG98+nUzG5Y8Hu+xbx6fTilyDZGyU5O4ZNl/ixxbmHic/Hf+5cTjZ9rRp2yNUhMm5XpNx+ba1QmdtdYhdJTaFlglJi33j+cxhGZsgD3MrNRe6MWmWZ1e5Ni8xOOufL/fLXLOxnH7Smw5SOXuH5vZ2xR+rVcTgn8t4VbFlYljudr1YuBPRcoiUnYK2JJln5eYLtn7eVQHrzW87SRLM7PRhGbpPSgcgJrompaujr7Ofma2jLunvbfzUvblJH8QdOX7/VmB/RB6ngN8WkL+H1M4YE+NeYwiNItfCV+0WkyKae51909KuI5I2ShgS5Z1ZOxy8t/8FoSaUina9eVsZrsQxhYPTuz+nND0/BJhDPbfgduA9dqTd4mSr/Mq4JJ2nFuw+byT5SjH+13sMx8Qt6X8ACqYj7svNrPrCb3o/5+ZreLu7xJ6m68Wk6mzmXQ7BWzpbWYlHr/n7u+X+wJxGNONhGC9GDgPuAHw/AlSzGxI6xzKIvk6G9vobd+Vuvz9TviU0HlvuRLStlXzn0II2DWEZvFfAv8bj80k1MJFupUCtvQ20xKPtwJuLZTQzLYkdH6aTui4VeyeedIkWjq6ne3uZ6Ulir2pV0g7VgZvEWrKdYTXWZSZnUTosf6mu/+9jOXojvc75zlCwF7bzAYXuo8dO5IV7ZPg7s+Z2XOE++J7EgJ2rrf7DUU6vol0GfUSl2pUcHrLMvhr4vGRbaT9OXA+obac/wVfrIxjE4+fLpgqDEPKjQcu649nd19MmBEMYAMz27ZQ2jjJy/mE8dI/KWc5KN/7XYrcmPb+wD5F0h1Aad99U+J2GzPbmZYhY2oOl4pQwJZq1JB7UO4m4zgm+h/x6Q5mlhqgzOx4wj1LCDW3fxYoY1r5kp2evpZyHDPbgqXvK5d1WFf0y8TjKWa2Sn4CM1uepYc5XVzOApTx/S7F9bTc+z7fzNZMuc4GhKliS81vMWHWtNxn9XKhcfUiXU1N4lKNPkg8PtfM/ki4D/tsmfI/FPgPMBQ4x8zGA38AZgBfIvQM/mZMuwg4PGVxjlwZR5rZjwljjOvjLGG3EGqqNcD34kQqf6HlHuuewH6EmmDOUDOrKeciIO7+oJldTqjZrkmYq/tXwCMxyWaEqVtzk5Tc5u63l+v6CeV4v9vk7vPM7BhCf4HlgafM7P+AhwmVkx0JC4Ikf2QV63z2qZlNJSw4slbcrdq1VEyvrmGb2Z1mdmelyyGt3EmYWxrCzFtPAWULJO7+OmGBkOlx106ElZeeINxjnUgItrOBPQvMzpW8F3tuLOPlMf/nCHNeE/M5mNBJ6QlC4D6QEKzvoWWM7wCWbkovl6MJU302E4ZK/ZQQwB4mTMGZC9a30vaCFx1Spve71GvdSPgR0kh4vWcDjxJq+acROgKeTEtv9YaUbJKuTjxuJEzxKlIRvb2GvebYsWPXpWPDg6SLuDuPPPIIl19+Oe5OU1MTyy+//KoLFixo3muvvbjtttsAePTRRz8olEdb6dydhoYGbrnlFh544AHcnblz5zJgwADGjBnDhAkT2H///YePHDny3laZx/Nvv/12pkyZwvTp06mpqWHs2LH/j/hvyd15/PHHufbaa3nhhReYM2cO/fv3Z7nllmPdddflm9/8JuPHj//6v/71LyZPngzAscce+1on37rUcgK8+uqr3HDDDfz73//mww8/ZPHixYwYMYKNN96YiRMnMn78+G/SejhX0dW72lmO5+ISn4cA36Blec2FwGvA3cBl7v5xGa51kZk9TOjlvT2hY9/nwGOE++SPE3ruQ/Gx5RBmPfscWAa4390L/psT6Wpduh62mX2X0JHlMHf/fQfOH0moqexJ6PAxm/Cf7sI4f3Fny/fS2LFj17377rs7m5VIT1S2gF1NzGwEYWgWhLXAjymSdhxhvW2A/WINXqQiuqyGbWabE37NdvT8FQjBeU3CL/8XCEF7L8JUi4e7+1XlKKuIZJ+Z7U/4fngDOM/d5xZI+tXE4+fbyPbguJ1FGW/
"text/plain": [
"<Figure size 525x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFICAYAAACSp82YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZyd4/3/8ddkzyTIRqkgJHxspWov/Qq11VpVWgSh1lpqqaWlihZFf9TWKkXUUluj1hZV1VqrtorwIUisQRYhmWyz/P64rjtzz8nZZubMnLln3s/HYx73uZfrOtc5OTmfc133tdQ0NTUhIiIiXVuvahdARERESlPAFhERyQAFbBERkQxQwBYREckABWwREZEMUMAWERHJAAVsERGRDFDAFhERyQAFbBERkQxQwBYREckABWwREZEMUMAWERHJAAVsERGRDFDAFhERyQAFbOlWzGyDPMcmmFlT/FuxGuXKKc9SZZTKMrNRqX/zq6tdHpFK6FPtAohUgpktB5wLHEMX/VzHHwu/BrYCVq9ycUQkY7rkF5tIG1wCHFrtQpRwC7AdMK3aBRGR7FGTuHQXvatdgDJkoYwi0kUpYIuIiGSAAraIiEgG6B62ZJqZnQ38POdYU3z4uLuPzZNmFHAKsDOwMvA5MBmYANzo7k25aVJp+wLjgX2ADYBhwGfAy8BdwA3uvignzQTg4NSh1VJlvNHdx+dc/yXgcOCbgMXnqAdmAs8BtwF/LlbOtjKz8cANcXcloAk4E9gN+DLhvXoW+J27P1BGfjsBhwBbAl8C5gNTgAeBK9x9Ron0qxPei7HAGoT3YgHwKfAMMMHdH2nNa4z59gL+CBwQDz0A7O3uC1ubl0hnUQ1bepq9gdeAHxICQH9geWAbQqB6wMzy/pA1szGEwHwNsAMhAPWN6bcHrgZeNjNra+HM7GDgHeAXhCC1UizjIGDVWP47gXvMrKPviRvwInAsMAroB4wAdgXuN7OrYuBbOqHZIDObCPwN+F4se39gCLAJcBbwlpntXvDJzX4CvAH8hOaA3xdYhvBvtz/wsJn9tg2v7WoUrCVjFLAl664GNgLuSx3bKP4dluf6K4Ea4CpCDXtb4AxgXjz/LeCk3ERxSNa/gXWARTH9rsBmwJ7AjUADsDbwmJmtlEp+VizP83H/o1QZz0o9x3aEWv5AYBZhmNq3gC0Igfqq+NwAuxd4fZV0G+EHw2Px+bck/ND5IJ7/IXBBbqIYxO8F9oqH7icE7c0JP2zOB+YAywJ3m9k38+RxSLyuD/A+cBqwYyzD94CbgcZ4+dGxJl8WM/sNodYOCtaSIWoSl0xz9+nAdDOblTr2UpEkC4Ht3P2p1LF/mtk/gKfj/sHARTnprgZWJDQJb+/uz+Wcv9fM7iIEqpWAS4Hvx/K8C7xrZnPjtYsKlPHcuK0HdnL3/+acn2hmf6P5x8k+wO+LvNb2WhG4Fjgy1fz+TKw5PwGMAU40s+vc/Y1UuuMJw9cAjnL33DI+ambXAU/G57jBzEa7+2IAM6uh+b34DPg/d38nlf4Z4A4zexa4Ih7bB3io1Asys/OAH8VdBWvJFNWwpae5MidYA+DuzwAvxN21083iZrYWsEfcPS9PsE7yuJ9Q0wbYx8y+XG6hzKyW0GQ8C7gvT7BOP8dncXflcvNvozeBY3Pvlbv7x8BRcbcv8IPkXKxdJy0Uf80TrJM83gZOj7urAN9JnV6N8D7MIfQJeIf8bk49LvlemNlPgZ/GXQVryRwFbOlpHixy7s247QUslzq+C6EZHaBUB6ck/16Ee9Blcfc6d9/U3YcD3y1x+fS47V9u/m30h9wOdAl3fxR4N+6m70NvQAjAUPq9+mvq8ZJmcXef6u4buvsQ4MdF0s8hdGKDEu+Fmf0IOC/uPoiCtWSQmsSlp3m/yLkFqcfp/xsbpR6/0Io+ZWuUe2GauzdC6LhFmMJ0NOHe+IbA1jQHxI7+wf1kifMvEDqTrWVmvWK50+/VJWZ2SZnPlfe9Sr0Xy8ZrRhP6EWxEeC8GxkuLvRc7A0ek9ucqWEsWKWBLT/NFmdfVpB6PaONzDW1tAjMbCZxMaIIvFPAb6ZzWsQ9LnP80bnsT3qNPqOB7ZWbrEJrXv0X+Ju9yh7WtFrcLCTXxfc3sVne/py0FFakWBWzpadoydjn9/2QzYHGZ6T4tfUkzM9uZMJZ7UOrwF4RhaK8SxmD/HbgbWK81ebdRfYnz6WFlSdN5+r36Ic0d+UqZn96JvcSvyclvFuG9mEQYC/4I8Dot369CHiaMvf8PIWj/zswed/fPiicT6ToUsEVKm5V6/IG7l6p5tlocNnYbIfgsJgyX+hPguZ2+zGxwpZ+/gGHAe0XOLx+3C1KBL/1ezSnRYz8vM/sKzcH6C+BsYKK7T825rhfNTeLFPAPs6e4LYi/xcwk9+f8fqQ5zIl2dArZIaZNSj7cAJha60Mw2J3Q2mwo86e7F7pmnjaO5o9sv3f3cfBeZWX/CBCKd4auEiWLylaMG2Djupq/Jfa9uLZS5mS1PqIVPBV5y9ySfI2n+bjrW3f9YIIuRlHdr4GV3T/onXAjsR7gPfqiZ/cnd/15GHiJVp17i0l00lr6kzdLje48uce3FwK8IteXce9DFyjgm9fj5gleFCUwGxMcd/YP7wCLndiZMVQotf8A8R3Mte/+4TnkhxxFqzxMIk88kyn0vxqUel/VexF7vR9B8a+Sa2LlPpMtTwJbuYkmv30o3Gccx0f+Ku9vH8bxLMbOTgW/E3ZcIM6PlK2O+8qXn1P5Wgfw3o3miEOj4YV3fNLOjcg+a2cpAMh3obELABZYExCvj7nDgptgqkJvHVoR7yhDuX1+bOl3Oe7ELqVniaMV74e5PANfF3dVpHu4l0qWpSVy6i49Sj883sz8CDe7+YoXyPwz4L2E6zfPMbBvCl/40Qk3zAJon/1gEHJFncY6kjMPjPNmPAHXuPpnQ2eynhN7pP4wTqfyZELxWJtRA9yNMVJJY1sxqOmIRkJTfmtmWhKbtOYSpQU8HVojnf+zun+SkuYCwWMjXCGO0X4rTgb5MmEt8e0JTeNJScLq7p//97qB5nu/z4zSvDxNmmRtFGKe+Fy178heryedzKqEn/grAcWZ2u7uX20FOpCp6dA3bzO41s3urXQ6piHsJc3lDaGp9DvhLpTJ39zcJC4RMjYd2BG4ndGiaSGiqriHUOPcsMBtauun4/FjG38X8XyKsikXM5xDCHNzPEAL3QYRg/SDNtdF+tGw+rrSrgLnxuf9G6PF9CSHI1ROmLL0+N1G8X7wj8I94aG3C1K5PEyZLOZnQWawBOMPdL89Jfy+h0xmE13xSfP6nCD8cvkN4jyYAydCsUfFHTlncfTZwQtztBVyXryVApCvp0QEbGD1mzJjdCfez9JfhP3d/8Zprrum90UYbUVtby4ABA1h11VVXnT9/ftNee+11cPIP/sQTT3xUKI9S17n7i//73/9GnXXWWWy11VaMGDGCvn37MmjQINZbbz2OOeYYnnrqqaHu/tcCZbzvwgsvZJ111mHgwIHU1taywQYb/F/q/Hk33HAD2223HSNGjKBPnz4MHDiQVVddlZ133plrrrkGd99lwoQJycIVnHDCCW9U8n284IILbqDZE4SZy64njMleADihOXxdd7+GAtx9prt/k1CLvYMwK9qC+PcGYQ70jdz9/ALpjyQs8vEIYVnRBsKPh9cJy2J+w90PoXle9b40LzZSFnf/E839E9YBftaa9CKdraapqanDMjezIwm/rA939z+0If1wQq1jT0KP0NmE2ZcuinM/t7d8r44ZM2bdBx4ouayvSI8wceJEfvKTnyS7+7n7bdUsj4g067AatpltSugx29b0XyJMjnACYRjL/wi1gL2AJ8zs0EqUU0REJAs6JGCb2VhCU9My7cjmdsK8wY8AI919E0LnntMJMyxdHacuFBER6fYq2kvczAYQAuqZtJy2sLX5jCV08JkL7B87iCQLAVxoZusTxmCeQcuxmCI9zqeffsqMGTNKX1jCOuvo969IV1axgG1mYwi9QlchdBA5Ezic5on3W2N
"text/plain": [
"<Figure size 525x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFICAYAAACSp82YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZwcVbn/8c8wCUwStpCwh5CQIQ9hFVGWC0hANoWALPGKBAlIkM2VRa5w2WQx4E9EiXABJQjIKsoSUBYRBQQRwg5PCJIACRGyISGBbPP745x2KkN3T89MdVfXzPf9evWrurvOqXq6CfP0OXXqnIaWlhZERESkvq2UdQAiIiLSPiVsERGRHFDCFhERyQElbBERkRxQwhYREckBJWwREZEcUMIWERHJASVsERGRHFDCFhERyQElbBERkRxQwhYREckBJWwREZEcUMIWERHJASVsERGRHOiVdQAitWJmW7v7823emwgcGV+u7+6zah5YQrEY642ZbQqcA4wEBgLzgJfdfY8MwxLp9pSwpdszszWA84ATqdN/82a2HvBjYGdgaMbhlGRmGwJPAGsl3l4XeCWbiER6jrr84yWSsp8AR2cdRDtuBPYApmcdSDu+SWuyvh24AvgIWJhZRCI9hBK29ASNWQdQgTzECLBF3C4BjnR3JWqRGtGgMxHpiH5x+66StUhtKWGLSEcU/mYszTQKkR6ooaWlJesYRKrCzM4Bzi6x+xF3H9l2lDjQBJwK7AtsCPwbeBmYCFzn7iX/hzGz3sBYYDSwNeFa73zgOcL13mvdfXGbOsnzt3Wdu49tU35dYBzwecDiOZYCc4CngJuB35aLs6PMbAjwRpki0919SCxbOO93gUnA5cAuhC70qcDp7v5g4tgd/s6KxPd5woDC7YB1gBnAncBFhO9mXix6lLtPbP8Ti9QntbBFWh1CGO18ArAJsAqwNrAbcC0wycyKjvsws2ZCkrkK2Iswcrp3rL8ncCXwnJlZZ4MzsyMJifOHhFuq1o8x9gMGx/hvA+40s6yviW8EPAbsDfQF1gA+TUjaQNe/MzPrZWZXAQ8CBxG+gyZgGPA94FnCjxqRbkEJW7qzK4FtgbsT720bH8cUKX850ABMILSwdwfOAD6M+79ASAQriLdk/RUYASyO9fcDtgcOBK4DlgGbAQ+b2fqJ6mfFeJ6Or99JxHhW4hx7EFr5fYC5hNvUvgDsSEjUE+K5AUaV+HydNTMRU7E4v1ikzncI92hfDOxKaEFf6O7T4ufpyndW8GNCbwOE0fXHATsBBwP3E3pIbuvcRxapPxolLt1WnARllpnNTbz3bJkqHwN7uPvjiff+bGZ/Av4WXx9JSEJJVwLrEbrP93T3p9rsv8vMbgfuIrSKLwW+EuN5E3jTzBbEsotLxHhe3C4F9nH3f7TZf4eZ/YHWHyejgf8r81krFruknwWoIM6ClQgJ+ozEe7cnnnf6O4txbA2cFF++DOzq7nMT9X9nZj8Fvl3BRxTJBbWwRVpd3iZZA+DuTwDPxJebJbvFzWw4cEB8eUGRxFM4xj2EViPAaDPboNKgzKwvoet7LnB3kWSdPMf8+HLDSo9fRVcUezOl7+woWm+FO75Nsi44BXipw1GL1CklbJFW95bZ91rcrkS4HlvwRUI3OsADFR5/JcI16Iq4+0J3/6y7DwAObad4YWrVVSo9fpXMcPe3S+xL4zvbP27fcve/FKvo7kuBq9sPVSQf1CUu0qpUgoEwm1dB8v+bbRPPn+nAmLJNKi2Y5O7LAcysH2EK02GE67zbEEZjbxSLZv1j/K0y+7r0ncUBdYXvr1y3PMAnekxE8koJW6TVBxWWa0g8H9jJc/XvaAUzGwScTOhOLpXwl5N9soZwbbqUrn5nA2n9jLPbqTOzk+cSqTtK2CKtOnPvcvL/oe0J9xtX4r2OnMTM9iUM2uqXePsDwm1oLxHuwX4Q+B2t04dmqdx32dXvbHnivYZiBRPK3sMtkidK2CJdkxzsNMPdU2/RxVugbiYk6yWECUFuArztBClmtmra56+Crn5n8wi3fDUSJkopp8M9GSL1SglbpGteTDzfEbijVEEz24EwcGoa8FiZQVltjaF1oNv57n5esUJmtgph8pF616XvzN2XmtnzhGvh25nZSoVr+0V8Kp2QRbJXD9e6RKqt1B/zNPwx8fz4dspeAvyI0Fpuew26XIzNiedPlywVJlBpis/r+cd4Gt/ZPXG7Lq0jxov5WoejE6lTStjSE3xceJJ2l3G8J7pwW9GeZvaDYuXM7GTCjF8QRjb/tUSMxeJLDqz6Qonjbw/8PPFW1rd1lZTSd3YFrWtwTzCzjYvUP5owe5pIt1DPv8JF0vJO4vmFZvZrYJm7T07p+McA/wBWBy4ws92AXxKmy9wAOJwwXSaEQVDHFlmcoxDjADP7H8L9yQvd/WXCYLMfEAZYnRAnUvktIZFvSJjK8zDCPNwFq5tZQ5qLgKSsS9+Zu79jZicS5ngfBDxtZuMJ85f3I8yKdlSNPotITfToFraZ3WVmd2Udh1TdXYRBSgDfJIyo/n1aB3f31wgLhEyLb+0N3AI8Qbg+ewgh2c4DDiwxs1fyOu6FMcYr4vGfBc6M+xoIieieePzfErp9exMmGSlMFLIyK3al15U0vrO48tYJhIF4AwhTxj5GmEf8aMKtZT+p4scQqakenbCBYc3NzaMIt6Do0U0f7j75qquuatx2223p27cvTU1NDB48ePCiRYtaDjrooCML/xgeffTRd0odo71y7j75+eefH3LWWWex8847M3DgQHr37k2/fv3YYostOPHEE3n88cf7u/t9JWK8e/z48YwYMYI+ffrQt29ftt56688l9l9w7bXXssceezBw4EB69epFnz59GDx4MPvuuy9XXXUV7v7FiRMnFhbD4Dvf+c6ULn53VRV/iIwgLI35APAvQvL9gHCt/jxgM3f/Q5ljXAFsBVwDvElojc8gzKO+Fa1zwIvkXlXXwzazbxAm+R/n7td0ov4AQsviQEK31zzCL+iL4/zOXY3vpebm5s0nTZrU1UOJdEft3eNc98zsUFpX7NJ62JJrVWthm9lnCSM8O1t/XeBJwjJ96wLPE371HwQ8GgeUiIiI9AhVSdhmNpJw68ZqXTjMLYR5kh8ABrn7ZwiDUU4nTJhwpZmN6GKoIiIiuZDqKHEzayIk1DNpXfquM8cZSRiQsgD4qrvPg/8sfDDezLYkTCZxRtyKSBHvvfces2e3N912cV/60pf+M+lIO2tfi0gNpJawzawZ+BNhtaBlhKQ9Dti4E4cbG7d3unuxvzZXEhL1l8ysj7sv6sQ5RLq9m2++mcsvv7yz1ZO3veX+erZI3qXZJT6IkKyfAHZw9wu6cKyd4vbREvv/Diwl3G/5mS6cR0REJBfS7BJ/G9jP3e9tt2QZZrYSrVMQvl6sjLsvMbMZhNb7cD45a5T0EMuWLWXu3Eqn5A5aWlpYtHBZ+wW7gQNGfYEDRhWdHK1dGw9pzn2r2t1vR70D0k2klrDdfSowNYVD9ac1rnJLEM4hJOzOrq0r3cDcuW8z79RhHaqzqGF1blzz+ipF1H1cfGndzrsi0iPV49SkfRPPPypTrnDdum+ZMkC437rEro79pRcREclIPSbsjvZVVn1GJhGRajCzIcAbJXa3ECaLehO4D7jU3cv1OtYlMxtLmPN9hrsPSrz/Z8LdQBe4+5nFa2fLzCYCRwKPuPvIbKOpz4S9IPG8qWQp6BO3C8uUAcDdtyj2fmx5b155aCIiVfMi8H7idS/CJcItCet6jzOzPdz9hSyCk+zVa8L+mLA84IAy5QrXrt+tekSSL6c/Qv8Bg0ruXuXDpXDZqyu8d/SxzfTpU4//O0gP8k13/3PbN+MUzdcRlgq93cxGxDkp8u5rhEuanZsooAequ79Q7r7czBzYGhhSrIyZ9SbMegYwpUahSU70HzCItdceUnL/gj4fAysm7EEbDWXVVet2CWnpwdx9jpkdSVjUZDhhZbOSC6Lkhbu/mXUMeVOvq3U9Gbc7ldi/PeHHxkesOLmDiEi34+5zCF3mELrIpQe
"text/plain": [
"<Figure size 525x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFICAYAAACSp82YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3debxd0/3/8dfNIBMhkgYVERI+EkOrran4Cj9TRaiS/qqCGKKKqg6GlqJaNfVHB1pftInqoKgWCS01tShVYgo+ESQIKRJzEpnu74+1trtzcqZ77z7Dvnk/H4/72Gfvvfben3Nzcj9nrb3W2i2tra2IiIhIc+vW6ABERESkMiVsERGRHFDCFhERyQElbBERkRxQwhYREckBJWwREZEcUMIWERHJASVsERGRHFDCFhERyQElbBERkRxQwhYREckBJWwREZEcUMIWERHJASVsERGRHOjR6ACkazCzrdz9iYJtk4HD4+p67j637oGlFIuxWZjZBGBSXD3Y3a9tUBzJ83b/5u57F9m/OnAWcACwPrAQeA0Y4+6z6hVne5jZ2YSYAXZw9wc7cI5ZwIbAve4+ukSZop8vM7sH2AX40N17t/faIgklbOkUM1sTOAc4nib9PJnZusCPgR2BjRocTm6ZWTfgdmCH1ObeQD9gTkOCagL6fEm9NOUfWMmVi4EjGx1EBb8DdgNmNzqQnNuTtmT9BHAmMBfo6+5LGhZV4+nzJXWhhC2d1b3RAVQhDzHmweap19929zsaFkk7uPvZwNk1vIQ+X1IX6nQmItXql3r9YsOiEFlFKWGLSLXSfy+WNiwKkVWUmsSlQwp63ibbkh7GRXvSmtkw4GRgb0IP43eBp4HJwNXu3lp4TOrYnsAEYBywFbA28DbwOHADMMndFxccM5m2XuoAG6ZivNrdJxSUXweYCPwfwOI1lgLzgIeBa4E/lYszK7E39knAgcAIYDmhVnsD8DN3f7fMsf2Bowj3nLckvI8WYD4wDbgRuKba+86pHtJpL5pZ8nqjzvYQN7PZwFDgSXffqkSZfwA7x9Vt3P0/RcpMBK6Iq6Pc/ZlqeonHjmMnAmOAjYEPCf/ml7j77SXimUw7Pl+p47YEvkW4770Obf8uV7j7X4odIwKqYUv9HAg8AxxH+IPYC/gYYbjLJGCqmRX9AmlmIwiJ+QpgD8IfuZ7x+N2By4HHLZVB2svMDickxB8Ao4H1Yoz9CInkQOB64CYzq/U9y82A6TGWTwKrA/2BT8Rtj5rZ0GIHmtlewCxCZ8Dki1EfQm/ujxMS0q+Af5rZGjV9F+0zNS63iF+cVmBm/YDtU5tGlzjPPnE5092fqebCZrYHMAP4DuHL4OrAQMLv729mdm4156nyWicSkvPhwAbAasC6wOeAP5vZFWUOl1WcErZ01OXA1sAtqW1bx5+ji5S/lFDLu4zwh3BX4HTgg7j/c8A3Cw+KNZ9/AiOBxfH4McC2wP7A1cAyQpK728zWSx1+Zoznkbj+WirGM1PX2I1Qy+9DqO2cE+PZnpCoL4vXBhhb4v1l6SzCl4Rb4vV3jNdM7hsPJ/z+VxC/sNwMDCCMj/4JsC/hfewPnAe8H4tvB3y3ynj2IfzO/je1bQxtv8tXqzxPOVPisoXwJazQLoQvaYnRhQXMbDVC6wiE30NFZrY14cvCGsASwu9sN2An4AzgHcLvaYMih1f1+UrpBfyU8Jn/EeF97glcEK8NMNHMvlhN7LLqUZO4dEicBGWumc1PbXuszCEfAru5+wOpbfeY2V3Av+L64cCFBcddTqiBvAvs7u4PF+y/2cxuIPyBXg+4BPhSjOcl4CUzS5LU4hIxnhOXS4G9ijS13mhmf6Xty8k4VkxetfAtd784tf6Amd1IuIWwLrC3mQ1299dTZb5HqLEBfLlI8+rNZnYdoam3B+F9fKdSIO7+NICZpSe+eTrjiVLuAhYAfQmtKL8r2J8k8SWExL2TmXV392WpMjsTEi9UmbAJXyR7Em45jHX3v6X23R8/W/cBgwoPbMfnK+1NYGd3fza17Q4zewz4Q1w/HLiuyvhlFaIattTLpQXJGoB4P/HRuLpZulnczDYF9our5xZJ1sk5phBq2gDjzOzj1QZlZn0JNZ/5wC3F7oumrvF2XF2/2vN30H8KknUSw1vA7+NqCysOs4Jwi+B14NFS90JjQpkeV2v9Pqrm7ouAO+NqsRp2su03cbkmoSabljSHzyck2bLMbAvgs3H11wXJOonLgVMqnasdvl+QrBN/JCRzgC0yvJ50IUrYUi+3ltn3XFx2I/whTuxDSEwAlcb8JufvRun7mytx9wXuvo27DwQOqlA8qWH2qvb8HXRbmX3PpV6vnd7h7nu5+zrANhXOX6/30V5Js/j6ZjYq2Whmg2lLYj+h7TbKrgXHfy4uby2oeZeyT+r170uWCp0NF1ZxvmoU/X8QOzI+H1fXLlZGRE3iUi+vlNm3KPU6/ZlM16AebUefso2rLZjm7svhow5OGxHuFW9G6Oy1E233MWv9Rbcjv6uPpN5Hb2AY4fdhhA5VOwKbxKItxY5voKmp13sQmv8h1K5bgLnu/pSZPUS4zzwauAg+GoEwMpavtjl8s9Trkk3Z7r7QzJ4k9JvorGr+bfV3WYrSB0Pq5b0qy6WTyEr3Das0oL0HmNkQwlCb/Sid8JdTn1apjvyuADCztQnDwQ4iJOli8dbrfbSLu88xs2mEL2p7EDpoQVtz+N2p5W6seB87qS0vBv5a5SWT3ujL4+2Gcv5b5TnLWVw49LCEZvsiJU1CCVvqpSNjl9Ofz21p60lbyRvtuYiZ7U0Y35yeyes9wjC06YROWn8H/szK941roUPjvM3s04Qm18GpzQuBZwm11f8Qkt25hF7ezWgKIWHvYmY941jxpOf33QXL/sCnCP8+ScK+x92r/cLz0e/ZzFoqjK/PYq70mo/fl65NCVua2fzU6znunsXwoRXEYWPXEpL1EsLQpz8Q+hu1FpRdPevrZyU2f/+JtmT9c+DXhIlIlhWUbdr3QUjY3yOMhd7ezF4jDHGDtkT9b8LwtNWB0bG5OrmfXW1zOIRhWBBaGwbS1umrGN1XloZTwpZm9lTq9faEGbqKMrPtCPc0ZwH3u3u5e4Vp42nr6PZDdz+nWCEz60VbE2oz2pe22cgmufuJZcoWnXSlSTxMaH5eh9AsnjwB62V3nwng7kvM7D7CeP7dgCcJw8GgfQl7eur1NpTo7BcfK7plO84rUhNNdx9Lcmd5Dc+dHmbz1QplLwLOJ9SWC+9Bl4txROr1IyVLhQlMesfXzfhFt6r3YWafJfXM5lKzyzVKbNVIelLvAfxPfH13QdG74nIn2ob+PebuL7fjcukvgEeVKTeGUAMvpZb/B0Q+ooQtnfVh8iLrptY4JvofcXV3Mys6M5eZfYu2OaYfI8yMVizGYvGlm0E/V2Q/ZrYtoYk50WzDoaC697EJcE3B5mZ8L8nwrm0IM4FB6YS9OnBEfN2e2nUy8cmf4uqBZjahsEwc0/+zCqcq9/kSyUxTfbuWXHot9fpHZvYbYJm7T8vo/EcTOkv1B841s10Ic2HPJsyNfQjwhVh2MXBMkc5DSYwDzew7hDHdC+IMXjcQpp5sAY6LE6n8iZAA1ydM6XkwK06L2b+KTkr1NoW2mcLGmNlfCHO0zyU0L+8FHMbKSWVN2sY1N4vbCf+WyTzb0JagE9OAtwgjApKWj3Yl7OhrhPvfawO/NrNdCWOy3yHchjk1xvABK3ZKTCv3+RLJzCpdwzazm82sI//Jpc3NhLm8IfzxexjI7IlD7v4cYR7pWXHTnoRZoR4kNGkeSEi2bwH7l5gNLd30+aMY4y/j+R8jzBlNPM8RhOT3ICFxH0ZI1rcCV8Zyq7FiE3TDxalij6Xt32J/wr/Dg8BNhIeurA48ROhYl6hHr/d2cff3gXtTm16IteF0meUFZea4e7lbGqWu9RqhdWY24d//MMKwsH8Rprldl/DQmXvKnKbk50skS6t0wgaGjxgxYixhuIV+OvDj7tOuuOKK7ltvvTV9+/ald+/eDB06dOjChQtbDzjggI8ePXj
"text/plain": [
"<Figure size 525x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"density = True\n",
"cumulative = True\n",
"histtype = 'step'\n",
"lw = 2\n",
"if perform_zscore:\n",
" bins = {\n",
" 'theta_energy': np.arange(0, .7, .03),\n",
" 'theta_peak': np.arange(0, .7, .03),\n",
" 'theta_freq': np.arange(4, 10, .5),\n",
" 'theta_half_width': np.arange(0, 15, .5)\n",
" }\n",
"else:\n",
" bins = {\n",
" 'theta_energy': np.arange(0, .008, .0003),\n",
" 'theta_peak': np.arange(0, .007, .0003),\n",
" 'theta_freq': np.arange(4, 12, .5),\n",
" 'theta_half_width': np.arange(0, 15, .5)\n",
" }\n",
"xlabel = {\n",
" 'theta_energy': 'Theta energy (dB)',\n",
" 'theta_peak': 'Peak PSD (dB/Hz)',\n",
" 'theta_freq': '(Hz)',\n",
" 'theta_half_width': '(Hz)',\n",
"}\n",
"# key = 'theta_energy'\n",
"# key = 'theta_peak'\n",
"results = {}\n",
"for key in bins:\n",
" results[key] = list()\n",
" fig = plt.figure(figsize=(3.5,2))\n",
" plt.suptitle(key)\n",
" legend_lines = []\n",
" for color, query, label in zip(colors, queries, labels):\n",
" values = lfp_results_hemisphere.query(query).loc[:,['entity_date_side', key]]\n",
" results[key].append(values.rename({key: label}, axis=1))\n",
" values[key].hist(\n",
" bins=bins[key], density=density, cumulative=cumulative, lw=lw, \n",
" histtype=histtype, color=color)\n",
" legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, lw=lw, label=label))\n",
" \n",
" plt.legend(\n",
" handles=legend_lines,\n",
" bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)\n",
" plt.tight_layout()\n",
" plt.grid(False)\n",
" plt.xlim(right=bins[key].max() - bins[key].max()*0.025)\n",
" despine()\n",
" plt.xlabel(xlabel[key])\n",
" figname = f'lfp-psd-histogram-{key}'\n",
" fig.savefig(\n",
" output_path / 'figures' / f'{figname}.png', \n",
" bbox_inches='tight', transparent=True)\n",
" fig.savefig(\n",
" output_path / 'figures' / f'{figname}.svg', \n",
" bbox_inches='tight', transparent=True)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFICAYAAAAoBEX4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xd873/8VdMEklMELmIa0Iin+K4pHVpiopSdepe1Z8SmihFkeo5FNVz0BYtp626q5aoUuE0p1RUSSjiEq26RfgEEVSFJK6DiMzM74/vd5mdnb32ZWbN7Jms9/PxmMeavdf3+13fvXeyP/O9rl6tra2IiIjk0Wr1roCIiEi9KAiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQiKiEhuKQhKt2NmW6c8P8XMWuPP8K6ul4isenrXuwIiCTNbC/ghcDz6tykiXUBfNNKd/Bw4st6VEJH8UBCU7qSh3El3nwhM7JKaiEguaExQRERyS0FQRERyq1dra2u96yCrIDMz4NvA7sAmhD+4FgOPAtOAG9x9eUx7FnBmSlH3uvv4mG4K8I34/HruvrDgeguAEcDP3P1kM9sbOAH4DNAIvAz8ATjf3d+JebYF/hPYDRga63c38CN3n9fBt6AsM9sIOBH4EjAS6Au8BtwPXOHuD6XkO4vwXr3j7mvHWbL/AewLbAx8DDwDTAUud/ePytShF/A14FBgO2AI8F7Mf0usR1NK3uSL47vAdOASYOd4/eeB09x9RkH6kcDJ8fVuDLwPPB6vcZOZXQEcA7zk7iNjnpOBC2IRR7r7NWVeyw7A7PjwRHe/JC2tSCG1BCVzZnYo8CQwGdgSGAD0AzYE9geuBf5mZut20vWvAG4D9iIEt/6AAd8H7jezNczsKOARYAKwASEIrR8f/93MtumMusX6fROYB5wCbA2sSXh/NgGOAB40syvNrE+FcnYGno7lfIrwPq8FfBb4BfCIma2TkncYIeDeCOxHeO19gcGEYHYB4GY2rsLL2Qh4ANiz4PqfJgTC5Fp7xXoeD4yO1xlE+ONjqpn9ntLzE34HLI+/H1ahHofH4zLghgppRT6hICiZMrPRwNWEL7oXgeMIX6rjCAEmaeFsC1wWf78CGAv8qaCosfHnqBqr8A1Ci8KBo4GdCF+QC+L5rYD/Ba4EXicE6nHAl4G/xDQDgV/WeN2qmNlE4NeEoPcioSW6C/C5WN/HY9JvAb8qU1R/4FZCAL2cEPA/R2j9/ium2Ro4p0Qd1gDuIbw3rYRgcyCwA/DvwMXAUkJgvNPMtixTj5MILcjz4+s4GDjX3RfEa3021nMAIUD9jBD8diG0aN8DDqEtiH0itvTviA93M7P1S1Ug/rFwSHx4m7u/Waa+IivQ7FDJ2mHA6kAzsJu7v1Rw7mEzuwn4K+EL+0AzGxK/7Baa2SdfXu7+OO0zhNAK3cXd343PPWhmjwFz4uO9CEFxR3d/I8loZn8BHga2Bz5vZmu7+9vtrMdK4pf4pfHhTGB/d3+/IMlDscv3t8DXgYlmNtXd72BlfQmzafcpOv+Qmd1OeK0DgEPNbLK7f1yQ5hxgC0Ir60B3v62o7DvM7LfAvYSu5N8QWpelrEYIemcUPPe/8fWuRugm7UMIgF909/sK0s2K/x7uI7TYS5kC7BOv83VCEC22F+Fzh9DLIFI1tQQla8lOLk20tUg+Eb+MzwQuIowndca/wR8UBMDkuk8TxroSZxcGwJimhbbWaC9gVMb1Op4QmJYDRxQFwKQOywmt53fiU98pU960UgHS3V8E7ooP1yR0swJgZmsTWpwAV5UIgEkZfye07gB2NLMdy9Tj8pTndyaMyUIYq72vOIG7P0sY00zzJ2BJ/D2tSzRpRS4C/lymLJGVKAhK1p6Nx7WAm81s8+IE7j7D3b/j7r8sDkQZaCF09ZXyasHvM1PSFNanMZMatdk7Hue6+0p/ICTixJ0H4sPPlxkbvLPMtV4o+H1gwe/jCYEY2gJlmtsLft89Jc2r7v7PlHP7FvyeOqmFMImnZIvb3QvH+MYW/3uKuwwl17m+qMUrUpG6QyVrvyVM1NiAMAlmfzObT/jCnQHMyLKLsYTFaTMagcKZkq9VkaZXNlUCM+tNGI8E2LpgdmUlA4B1gVKBZkGZfIXvQeH/87EFv08Lk3irsmnK86+UybNtPL7j7s+lJXL3j2N39W4pSa4hzKSFMK5c2PV6MGF8FdQVKu2glqBkyt3fAvYgjK0lNiVMVrkZWGxmM83ssDhFP2vvVZMoWZ7RhQbR/v9vg1KeTwv2ECa8JArf5yHFCTtYh3dTnocQvKGtO7Oc1B4Bd3+MMM4LYYyz8PUkXaFPdmAcWXJMLUHJXBznGRfHkQ4izLxMZhg2AF+IP5PMbF93/zDDy3d1cKtW4f+1O4DTa8j7fOUk7arHAcBLaQmLvJPyfLkWbd94rCb4V2oZX0NY9jGSMKt1lpmNIMwyBbUCpZ0UBKXTuPtswgLm78VF3V8gzPQ7gDDFf3fCAuof1a2SXadw2n5DHVsthfVY1Mn1WExYnznEzHq5e7lAV6mFej1hok4fwgL/WfHYi/CHz/Udr67kkbpDJVNm1t/Mti1eW+buC939Bnc/lPCXfEs8tU+XV7IO4s4tSYvuM3GMMJWZHWtmx5vZ3pXS1mhOwe9pyx6SOowxszNi1/Vm7bhWEmAbgdT8cSnFtmnnAdx9EWFnGgiL+yGMOQPc4e6vt6N+IgqCkh0z60sY/3mM9GnzyRhPMjGlX8GplhLJVyXJYvx1aFvcvZK44cClhDV2F2c8fjmTti7joyrsSvMD4MeExfSVdo4ppXDzg3I7vuxNdWOVU+JxhJntWVAndYVKuykISmbidPZk2v7OZnZQqXRmtith9ijA3wpOfVSQJuvlCd3BLwmbCABcaGZbFScwswGEoJP837woywrEjQmSJQebAxeXmqBkZgfTFrgWEiY11epOYG78/ZS4v2fxdTYk7FBTjem0TaC5hPAevcWKwVakJhoTlKydTfjLvjdwY9x5ZDphjd5gYFfCxtoAH9C2QTKsuGzh3Ji3ObYcezx3f87Mvg/8lPBezDazSwkTZZYSllD8B21dh4/QtsNMlpJNwzcizNrd1swuI2w1N4zQ3TiREGRagWPbM3nJ3VvN7BjCus3+wF/N7CLC611G2DXoe6y4W0zquKG7Lzez6wmbLCTv0Y3lNgkXqSTXLUEzu9XMbq13PVYlMWB9g/Cl3ptwp/g/EJZMTCd86TUSJmh8xd29IPuttLWUTiS0Ev/YNTXvGu5+PuE9WE4IDCcT1k/OInQhJ1/u9wNf7ozF3+6+GPg88ER8akdCl+LDhM/gKMJn9yEwyd1v6cC1ZhHW9n1IeL2nEoLiA4Q/gIYS9pBNxksrBbTiRfdT2ls3Ech5EARGjR49el/CX5/6yejH3a+fMWNGvyOPPJItttiCxsZGevfuzaBBgxg7diwnnXQSs2fPXidu+VWY77Ff/epXDWPHjmXAgAH069ePjTfeeOMPP/ywFWg98MADk9soMWvWrNcK826wwQYjADbZZBNLq9f48eOTHVtIS3Peeed98iX729/+9p5Oen/OnzFjRu9JkyZhZgwcOJDevXszZMgQdt11V372s5/xzDPP7BKD1Qp5TzjhhDOT+k2dOvWhtGtUSufuL86dO3eb888/ny984QsMGzaMPn360L9/f8aMGcOkSZOYMWNGf3efknKNqrn7VEIr9zLCTjZLCTvEzAD2c/fjCbM+ofzaR9z9Kdq2v3vW3R+ppS4ixTr1foKxK+QK4Gh3/3U78g8mDM7vT7gNz1uEvyDPd/eHy+WtsvynR48evcX06dMrJxaRYlnuqPMuYXu3P7n7fmXSrU0Yo1wdON3df5JVHSSfOq0laGbbs+J4T6351yWsMTuJsPPEk4S/QA8kLJQ9Mot6ikjnMLPdzeyPZvY/ZrZxmXTb0ba/6RNp6aJDCQFwOWGLPpEO6ZQgaGbjCdPBB1ZIWs5Uwi7
"text/plain": [
"<Figure size 480x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFICAYAAAAoBEX4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhcVZ3/8XfoBJIQlpAQEhMJJCFfCLJkZDEEhqCAKKtsowISGBBcBnVk0YFnABECMj8HERgUl7ALKD9hCKgQAhKWsBiQzW+AEFAQgbAGEpJ0ev4459JFUWv37a7qOp/X8+S5VX3PPXW6q9OfOveec26/jo4OREREUrRaoxsgIiLSKApBERFJlkJQRESSpRAUEZFkKQRFRCRZCkEREUmWQlBERJKlEBQRkWQpBEVEJFkKQRERSZZCUEREkqUQFBGRZCkERUQkWQpBERFJVv9GN0D6LjPb0t3/XOLrM4HD49NR7v5Srzasm8xsGjAnPv2uu5/doHYsAsYC7u6bltg/ADgROATYEGgH/g4c4e739l5La2dm04FfxqdfcPdfdaGOO4CdgefcfaMyZVryd1PypxCUupnZOsD3gK+h36FGugI4uOhrawMvNqAtTUG/m1Iv/ZJIV/wQOLLRjUiZmW1KZwA+B3wHeBYY6u7PNaxhjaffTamLQlC6oq3STnefDkzvlZaka/OCx2d25bRiI7j7TGBmD75Exd9NkWIaGCPSN61Z8PjZhrVCpI9TCIr0TYX/d1c2rBUifZxOhybOzAz4KvApYGPCH9dXgYeA64Gr3H1lLHsacGrR8R3x4Z3uPi1+bSZlRuAVjHj8f+5+vJntCXwd+DgwBHge+A3wA3d/Mx6zNfBtYBdg/di+24Ez3H1BLj+ICuIozK8AXwA2I/y/eR64ETiv0ghDMxsIHAZ8FpgMDAcGAK8DjwE3AZe4+zs1tuUOwsjIQnPC2wjALu5+Ry11VXiNPwI7AW8Aw9x9VYkylxG+L4CD3P3XJcrsBvwhPv2su99Sy+hQM1sXOAY4CBgP9AMeAS5y92vKtPk0avjdLHHcRsAJwB7AaOAt4AnCKdtL3b2j1HHSOtQTTJiZfRH4M3Ac4RrTYGAgMAbYF7gUeMDMNuih17+YEAJ7EMJtEGDAfwB3mdmaZnYUcD9wKOGP1OrAR+LzB81sq55oW4FRwIPAj4BPAOsQTkVuBpwEPGJmk0sdaGb/BDwF/BTYjxD+a8bvYQPCB4//Buab2eie/TbqclPcrgtsU6bMpwoeTytT5rNxu4TwoaWq+IHnCeBswgejdQk/838GfmVml5Pf360DgCcJHwLHAWsQfg93JgT1LDNTR6HFKQQTZWYTgF8Q/iA/S+jp7AhMIQRMNs9sa+Ci+PhiQm/mfwuqmhz/HVVnEw4nfNp34GhgKqFnsSju3wL4NfAT4B+EoJ5C+MP6+1hmLUI49aTjgC2Buwjz8XYg/HweiftHAFeYWb/Cg8xsGHAr4QNFO/AzQhBOAfYETib0aAE2Af6rxvYcRfh5F/Z6jqbzfXiw9m+trJsKHu9WvNPMJhE+iGSmlaknC8E/uPt71V40fhC4k/DBo4PQG9uD8DP/BmEO5KGE35ViXfndvIDQy7wwvs4uhPcl65V/Bvj3au2Wvk2fctJ1COGTbzvhFFrhsPr7zOxa4A7CH6DPmdnweNrvJTN7LSvo7g938fWHE3qhO7n7W/Fr95jZfMJpQgh/mBYB27v7y9mBZvZ74D5gW+CfzWxdd3+ji+2oxfnANwtOjd1rZtcBDxACchKhx/RAwTHfAtaLj4939/OK6rzZzC4FHif0dPYzs/7Zqedy3P1peL/HlHm6G+9Dqdd4wswWEnpHuwFnFhXZNW5XEE7tTjKz9d39layAmY0HJsanN9b40j8gzHME+LK7/6xg371mdg3hw8gmJdrcld/N94BPuvs9BV+7w8xup/ND4OGxXdKi1BNM18i4XUKJydXuvoLQ2zif8Ae9J35XTikIwOx1HyecosqcXhiAscwqOj/x9yNcN+opLwInFF8bcvflhN5d5mNFx42Ox75E6Gl8iLu/QPigAeE09LAc2puXWXE7xczWLNqXheDlhB5bPz58nTLrBbYX1FVWvA6YzXu8rSgAAXD3fxDOHuTlgqIAzF7nPuBP8emmOiXa2hSC6fpL3K4DXGdmmxUXcPfb3P0b7v6j4iDKwSo6lyYr9kLB49llyhS2Z0guLSptdgy8Up4qeLxe4Q53P8LdRwOj4weKcgoH1azRxTb2hOyU6OoUBFwMhOz5r4FsYNIuRcd/Jm7vdfdXqW53Os9MXVWukLvPIb8pITdX2Je9t6sR/o9Ii9InnHRdRhgVN5owCGbfeArsVuA2wqfxnjzF+Kq7Lymzr/D60d9rKNOvTJk8/K3CvmUFj0v+X8pGVsYRpmMJpxgnEq55TonbTDN9KL2DcJZgCOGUaBYY2xFOWa4knJq8mzCYaVp2oJkNKnhe66nQwrVRq53GfIAwkrm7uvXeSmtopv900ovc/XXCaa37Cr48jnC66TrgVTObbWaHFA/6yMnbtRSqdo2sF9TUTkoEsZkNNrPjzexPwLuE3sXvgR8DXyYE4IemHzSD2Pu9NT4tHByTnQp9IH6IyXrzk8xs/fh4F8JIX6g9BAtHIC+uUvYfNdZZTZffW2kdCsGEuftf3H0KYej/uYRBGpk24JOERZpvjZ/u89TocKtVl+aJmdnGhBGk5xJGKPYHlhMG/VxHWOtzCmFUY7PKTolubmaj4uMsBOcUbaGz95ddD1zg7l7jaxX+nKuFTqXTy/XQHEBRCAq4+zx3P9HdP0YYnn4IcDWwNBb5FHB8o9rXR10NTCh4vAMwxN23cPeD3f2cOACjeNBJM5lFZ1DsamaDCR+YIIZfHNyTXRecFrfZ9cBae4HwwdPe65ctFaxXZb9IzXSuO1GxZ2fAijgiE3h/qPlVwFVxEviDhA9LewFnNKKtfY2ZbQNsH5/OcfcvVii+YS80qUvc/R9m9iBhKspuhNGuAwg92rsLit5OuM75yXh3i3Hx6/WEYOFZiG2pPN9x6wr7ROqinmCCzGx1wnWX+cD/lCvn7vPp/IQ+sGBXU17HaiITCh4/VK5QXLKrcOJ3M34ozU6J7krnqND73H1pQZlsNZhN6byN0avAh6YfVPB7OiepH1HuOrSZbUHlENTvptRFIZigOOghW9NxRzM7oFQ5M9uZMHoUPjgR/L2CMj05PaGvKpwSsGupeWZmNpJwbXD1gi830xSJTBaCo4Aj4uPiqS1z6Dxt+m9xe7O7t9f6IjFUs+uj21K0DiiAma1NWOWoEv1uSl2a8ZOn9I7TCct39SesyXgZ4RrQC4RJ2zsT1lSEMLLx3IJjC6/fnBWPbY89R4G5hJ/RKEKv5XYzu4Bw89v1CD/bfyWsmlOo6eajufufzOxFwjJpY+KXby8q86qZPUpYPSc7Y1DPqdDMqYSl5cYDp8a1V39GGA26BWGt1gmEHmO5a6n63ZS6JN0TNLMbzawr/1n7vPhH4XDCfKj+hNNYvyFMmZgFnEiYI/YasH/RKL8bCSuBQPjk/wDw295pefNz92WEdVCzU4Y7AdcQfrY3E/6YDycMKDmp4NDCG+U2k8IVX5bywWk1mcJgfI/O9V1rFu+k8c90rsu6N3BDfL1LCAF4M2GlmnL0uyl1SToEgfETJkzYm3AqJ7l/7n7lbbfdNvDII49k0qRJDBkyhP79+zN06FAmT57MN7/5TebNm7eeu/+u6Lj5P/3pT9smT57M4MGDGThwIBtuuOGGS5cu7QA6Pve5z2W3UWLu3Ll/Lzx29OjRYwE23nhjK9euadOm7VnwHpUsM2PGjOx2PFx22WVzypXryr9YHwDf/va3Z3SlnLvfdssttww66KCDGDNmDAMGDGDAgAGMGDGCqVOncsYZZ/DII49MfOihh84ZNCjMPtlhhx1+WuvPqie//+J/F1100dHZa02ZMmVQXAy7gw8qDME5FRZCqMjdXyQMKjqKMPhmMaHnN5+wmPnedIZcqeMfBvYhrP25hBDaK+LIVpEP6dfRUfy
"text/plain": [
"<Figure size 480x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFICAYAAAAoBEX4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dedxc4/3/8VdWScSSxB5CFj5ELamtQb9SLVVUKPqtIkKtLVq19suX0KLor9WW1lKEoEWtbbRq+0ZjiVZR6yeCxFZLYo0kEsn9++O6jnsymfWeM/fMPef9fDzux5kz5zrXXDOZnM9c6+nW1taGiIhIFnVvdAFEREQaRUFQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQREQyS0FQmpKZbVrk+Ylm1hb/1ujscolIa+nZ6AKI5DKzlYCzgO+h76eI1JkuMtJsfg4c0uhCiEg2KAhKs+lR6qC7jwfGd0pJRKTlqU9QREQyS0FQREQyq1tbW1ujyyAtyswM+C7wZWAo4UfXbOAx4Bbgenf/NKadAJxRJKsp7j4mppsIHBSfX9Pd38x5vZnAusD/c/cTzGw34GhgC6A/8ApwM3C+u38Qz9kcOB74ErBqLN99wI/dfXqNH8EyzGwMcH/cHQ08A5wC7AMMARYAjwNXED6fuvwH7czPysyGAocBY4BhwEDC+3wHeASY6O53553zOeCfwHLAJ8Dm7v58gbx3BO4BugHTgc+7+8dVfyCSWQqCUhdm9m3gKqB3iWRPALu4+1tpB0HChfyIIvk9RQhA+wG/AXoVSPMR8EV3f7JE+auWFwTHAucBGxZJfhuwn7svSLMMsRwz6YTPysx+RBjtW278wW/d/bt5554InB93pwL/lfujwMxWjuVbG1gEjHb3x8q8jshS1BwqqTOzEcCVhAD4MnAUsD3hYnoA8HBMujnhwgpwCTAK+FNOVqPi36FVFuEgwkXdCTWQ7YADgZnx+CbAH4FLgbeAY2PZdgXuimlWAH5Z5etW61JCAHwc2B/4AmHQz3Px+J7A1XUuQ90+KzM7GDiHEABfA04Gdo7n/zdwLbAkJj/KzL6al8X/IwQ/CN+fI/OO/4YQAAFOUwCUjlBNUFJnZmcAE4DFwHB3n5V3vBfwf8C2QBuwmrvPjscmEmt67t6tQN6fHad4TRDg34TayYc5xzcGns7Jbiawjbu/nZOmO6GJbqtYtoHu/n7Fb76MvJogwJ3AXu6+MCfN8oQmvi/Ep3Z099xz0ijHTOr4WZlZN0KT6trA+4RmypcLlONo4Ndx9wp3PzTv+DDgSUJt9UNgpLu/bmbfAn4fk90PfMXdlyBSJdUEpR6SlVzmAm/kH3T3RYSmz18Bx1Gf7+FpuRf1+LrP0F7LAjgz96Ie0yyhvTbaDRheh7Il3gMOzA2AsQwfA+NoryXl14DSVo/Pal3gXeAD4KpCATC6Nufx4PyD7v4ScELcXRG4MK4UdHF87l1gnAKgdJSCoNRDMoBhJeAmM9soP4G73+Pu33f3X+ZfXFOwhKVrW7lez3l8b5E0ueXpn0qJCvuDu79b6IC7v0B7U+AuZlavOb11+azcfaa7b+buK9MexAr5AJgfHy9XKIG7Xwr8Ne7uEx8PjPuHuftrJfIXKUmT5aUergFOJPyyHwuMNbOXgLsJzXz3pNnEWMBsd59b5NgnOY//U0GaZZpkU/RgmeP/Av6LUANai9C8mLa6f1ZJLc3MViSMDh0ObETo790e6BuTlvpR/h1C8+wAYLP43BXufkuJc0TKUhCU1Ln7e2b2FcLo0KRfaxhhAMYRwGIzm0IYPFOPaQAfVVjOT1N+3Wot01Sc552cx2tQnyBY188qtgL8EPgaBZo7CX2Jlbz+G2Z2Gu3NoIuAkzpSJpFcag6VunD35919NCEIXkCYD5foAexI6A+628z6FsiiFo0ObpUqV87cJeQWFk1V3zJ0WBwd+m/C6N4kAL5LqAFfSlgjdghQdl5fHISzX85TvVAQlBSoJih15e7TgGnASXFAw47A7oTh/30JE+lPAH7csEI2zsAyx1fNeVysObIpmdkmwGWEa8xHhNHCt7j7zLx03WlvDi3lBELTKYR+xJWAE83sDnd/KKViSwapJiipM7O+ZrZ5HGb/GXd/092vd/dvE+ajJSP6du/0QjaHzcsc3ypu/8PSA1C6giNo/5F9tLv/PD8ARmtT5joUV485K+4+TfjufBLPu9rM+qVSYskkBUFJlZn1BuYQJoD/tlg6d3+c9tpNn5xDWRrqvp+ZFbxrhpmNBLaJu7fVa/m0OhqR87jUJPYDch4v0zIV55ReQxg5ugT4Tpy+kbQcjCA0t4t0iIKgpCrOeftb3N3ezPYulM7MdqC9n+gfOYc+yUlTz+kJzcCAs5d5MoyivIYw2nIRcFEnlysNs3Mef61QAjPbFTg956lCUyROJ4wiBfiluz8aH59P6G+EsNrMV2ooq2SY+gSlHs4EdiN8v/5gZtcAkwnzzgYBOxAW1gaYx9K/5HP7vs6J5y6ONcdWdHJsNr6M0OS5KWFB7WHx+E/d/dlGFa4GNxKWgoPw77gm4cfRh8B6hPl+e7H0tIqVcjMws62BH8Xdl4HTkmPuvsjMvkNYsaYHcKWZbZIs9i1SqUzXBM3sDjO7o9HlaDUxYB1EuFNAT8IowJsJF6zJhFF9/QkjBb/h7p5z+h2E5dYAjiHUEm/rnJJ3uquBNwl9oncQPp/LaA+AZ7n76UXObWrufgfhvUAYyflDwiT3h4DrgW8QAuBE4PaYbr2kfy+OGL6G9hGyh7v7vLzX+CdwYdxdh7ACkUhVMh0EgeEjRoz4OmGukv5S/HP36+65554+hxxyCCNHjqR///707NmTAQMGMGrUKH7wgx8wbdq0ge7+17zzHr/ssst6jBo1in79+tGnTx+GDBkyZP78+W1A21577ZWsG8rUqVP/k3vu4MGD1wUYOnSoFSvXmDFjdsv59y+Y5txzz70qSXDNNdfcn+bnEvMD4Pjjjz/o4YcfXmPcuHEMHjyY3r17M2TIEMaOHcutt95KDIB1+fdJ6bMqyd2PICyUfTehn3gxYSm95wkB7ovufjDtS6/1ItQOAc4lNBdDWHbtniIvczrwYnw8zsz2LFcukVx1XUDbzI4g3B3gMHf/XQfOH0RoAhlLGEX2HmGO0fnu/kgK5XtmxIgRIydPnlxrViIVmTZtGuPGjQPg+OOP5/DDD29wiWpSz9V0RDpF3WqCZrYVNYzaMrPVCfPLfgCsTugEbyP8UpxqZoekUU4REcmuugTBeLuYuwj3GeuoGwhrDN4NrO3uWxLWTzyF0E9wSaGFmUVERCqV6uhQM+tDCFKnsfSST9XmM4YwgnAu8G13fw8+W4j3vDh59gDgVJaeZySSuvfff5///Kf2BVuGD6/trkzPPfdc+URlrLLKKqy66qrlE4pkRGpBMN5N/D7CKK3FhEB4GO037qzG+Li9PbnZap5LCMFvTzPr6+7zC6QRScV9993Hj370o/IJy7j33mJ3I6rMnnvWPubj6KOP5phjjqk5H5FWkWZz6NqEAPgI4Q7Uy0wCrsLouJ1a5PijhIV/lwe2rOF1REQkw1IbHRprghu4+505z80k1AQrHh0aF9T9hFBL3dnd7y6SLsn7UHe/ooNlrvvo0EVzXuXDqRPrlr9Iowwa+78aHSpdXmrNoe4+A5iRQlYDaC/XOyXSzSEEwVVSeM26+fS915hz64RGF0MkdYPG/m+jiyBSs2ZcNi13RfgFJdIl/YBlV5A3s2eKHKptpIKIiHRpzbhizOLySZZSv9n+IiJNyMyOMLM2Mzu0yvO+F88b34HX7BnPbTOzsqPyzezQmLapb3LdjDXBuTmP+xRN1X4jznkl0gDg7hsXej7WEEdWXrQUdO9B/1F7dOpLikjr6OhCJGa2DXBe+iXq2po1CH5CuK3KoBLpkr7ALnWz0e69+7HWMX9sdDFEpAuKc6hvocqFSMzsy4RF7JevQ7G6tKYLgu6
"text/plain": [
"<Figure size 480x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFICAYAAAAoBEX4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df7xd053/8Vdyg4QEESKakEgiH9KqapVm8JUpozqUKjolfoShaNW0MyjlWzEtxo+ZL0UboyVoU6RNUb9aPyoahNLE73ySiIimDRKhgojce79/rHXkODm/z7733HvX+/l45LHPOXvttdc55+a+7957rbV7tbe3IyIikqLezW6AiIhIsygERUQkWQpBERFJlkJQRESSpRAUEZFkKQRFRCRZCkEREUmWQlBERJKlEBQRkWQpBEVEJFkKQRERSZZCUEREkqUQFBGRZCkERUQkWQpB6XRm9skSr08xs/b4b0hnt6urKPX5pMrM1jez7Yu8PiLv52VyM9om3Z9CUDqNmW1iZpcDf252W7oiMxttZncDP2p2W7oKM9sHeBr4WrPbIj1Tn2Y3QJLyP8BxzW5EF/Y7YCQwo9kN6QrMbGvg3ma3Q3o2haB0ppZyK919IjCxU1rSNZX9fBKkz0M6nE6HiohIshSCIiKSrF7t7e3NboN0Q2ZmwDeAvYFtCX9QLQOeBKYDU919TSw7CTi3RFUz3H18LDcFOCa+vpW7L83b3yJgOPDf7n6ame0PnAJ8BugPLAZ+DVzs7m/FbT4F/Afwj8AWsX0PAD9w93kNfgQlmdnehGuf/wB8DFgNLAVmAte7+4MF5R8E9ipR3XnuPimWW0T4DC4HLgSuAPYDegGLgAvdfWpevb2ArwJHALsAmwNvAy8AtwGT3X1lifeQ+8XwHXe/zMy+AhwPfBoYCLwK/AG4zN1nV/g8vgScCOwKbBK3vRe4yN3nmdkqYIOC91ruF9M/uvuDZjYCeCm+djVwMnAU4ZT6TsCGwBLgntjOBeXaKWnSkaDUzMyOIPTYOxX4OOGXTV9gGHAQcD3wJzPbsoP2Pxm4gxAAWwD9AAO+B/zRzDYys+OBx4EjgaHA+oRAOhJ4wsx26qC2XQncRwieEXG//YHRhF/OfzCzG82skevxGwMPAYcBA2L9nyD8ws+1YzDwR+Am4EDCe18fGATsAVwCuJmNq7CvFjP7BeEPjC8CW8Z6tgaOJnyWJxTb0Mx6m9n/ArcD+xO+q9y2xwFPmdkhtb75EgYQgvV6wh89mxF+JkcB3wRmm9kXMtqX9CAKQamJmY0GriX8MnuJ8Nf3HsA4QsA8Got+CvhxfDwZ2Bn4bV5VO8d/x9fYhGMIRxUOnADsTvjrf1FcvyPwK8KRwauEoB4H/DOh9yWEX5iX17jfiszsKMIvXAg9PA8HdiMc5Z1KOFqF8DmdnLfp8YTP4m/x+ZOs/XyKjX87BhgD/CzWfSBwhbvPiO3YiHCUtjvQDvwcOJhwJPZFwhHkKkIw/t7MPl7mbZ1GCPQXgK8Dn4t13BzX9wauMLNtimx7OeE7Angxb/sDCaHaF/gl4Wep0M6E4My5mrWfyRNFyh9BOCvxNOHz3J1wFPyHuL4/cIOZDSjzXiVB6h0qtZpAOHXVSjgt9XLeullmdgvwIOFU4MFmtnk8rbnUzN7IFXT3OXXuf3PCL7o93f3v8bVHzGw28Gx8vh8hFHdz99dyG5rZ74BZwGeB/2Nmm7r7m3W2o5jc8I/ngX3dfXXeuofM7Dex7QMJgXAFQO40nZnlyq+s8Pn0Jpxuzv8DIv8PjPOBscAa4GB3v6Ng+3vM7AZCUPcnhOnnSuxrCHA/cIC7ryqoYwVwEuHn4XDgotxKM9uZcLocwhH5Pu7+dn57zew/gEuL7dTd55hZ/neztIqfmenA4fmfu5n9CriTENyDCT8b0yrUIwnRkaDUKjeTy0rgr4Ur3f0DwvW/HwHfoWN+xs7JC8Dcfp8jHK3knJcfgLFMG2vDohfhVFmWcp/NywUBmNv/X4DvAxcTjp56NbCvnxR70cw2Ze3R1zVFAjDXlidiOwB2M7Pdyuzr1IIAzLk673Hh6eVvEb77duCYggDMteG/CQGbhfeBrxd+7u7ezkc/q09ktD/pIXQkKLWaG5ebANPM7Cx3zw8f3P0+wnWxjtDG2lNchZYAO8THpX655gdj/6waFc0Ftgf2M7PzgcuLBPGVGexnDcVPCQKMJ1yjhcoDze8CJsXHewOPFSmzxN2fL7H9i3mPPzzNaGa9WXsq81F3n0tpV8d9N+pxd19eYt38vMebZbAv6UEUglKrG4DTCZ1NDgIOMrOFhF+49wH3ZXyKsdCyUj0aCUcDOX+rokwjR2LFXAocQPh/9T3gTDP7M+FzuReYWewIsQ7LShyZQbhmljM9dOKtysgSry8qs03+95D/u2QrwqlHKB3WOcWCtx5/KbMu/7PS7zz5CJ0OlZq4+wpgH8K1tZyRhM4q04BlZna/mU1o8HRfKeucVismNzyjM7n7w4QOKLlfyL0JQxPOJByZLjezm8xsjwZ39fcy6zavs86BJV4v9QdH7lRjTv53nd8reFmF/b5WYX21qvq5IPs/fKSbUwhKzdx9rruPI3SmuAR4Lm91C/B5Qo/Ee82sX8a77/Rwq0W8BjeS0APyWvKGLRBOv/4LYRjHBQ3sptwYuvwjnS+ztkdlpX+nNdCeQvm9PSv9jslqoLIGPEtddGpA6ubujxFOZ50Rb330ecLpwC8Txu7tTfjl+oOmNbIJYueg38Z/xNsA7UP4XD5POBo5y8zucveZGe/+jbzHrzfQC7cR+Ud/W1QoW++Rq0gmFIJSk3hkZ8AHsUcmAHEYxFRgauwe/wThKOAAEglBM9uM0DHmJXf/8Jpk7BgyF7jSzE5l7RjFLxFmkcnSs3mPPwc8Uqa9YwgD7hcROpbML1W2RgsJp2w3JswwU85nMtqnSF10OlSqZmbrA8uB2ZToog8Qp9HKhUDfvFVtHde65orX+ZYDD/PRgfCF7sp73LdgXRafz/2sPWV8vJmtV6bsOcAPCaeuK80cU7U4FCX3Pj9nZuWGohxTZl2P/XmRrkMhKFWLPRt/H5/uUWrKKzPbi9B7FOBPeavezyuT9fCEZnuctZ08Tokz6xQzIe/xnwrW5T6fuj+bvCNyCMNFio5HNLPD8tqylOwHkF9OuE7XC7jOzDYsLBBn2PlKmTrye/L2tJ8X6SJ0OlRqdR5hDFgf4KY488idhA4ggwjTeOVmCnmX0HEmJ3/YwgVx29ZKEzB3B+6+2sx+QJgFZiDwuJldQbhm+gZhirJDWXuH9LmsnXos52+E06k7xblP5wAr3P1FapObNHxrQq/dT5nZjwlTzQ0mdNqZyNrB7Ce5+3s17qMsd59lZlcTZpTZE3jSzC4lnK4dSOggVHgUWNi5ZRnwAbAecLiZ3Qu8Ccx39zcQyUDSR4JmdruZ3d7sdnQnMbCOIYy96kOYKuzXhCETdwJnEP5qfwP4irt73ua3E6ZbgzCjyJ+AWzun5Z3iKsJMORB+0X+f8Jk8SviMDiccGT0P7Bc70OSbHpd9gGsIn8/3a22Euy8D/g/wVHxpN8LE0rMI38HxcR/vAce6+2217qNK3wJuiY+3B34a23A3IYTf4aPvL//ID3dvje2F8EfEPXH7gzqovZKg1I8ER40ePXos6l5dE3fnlVdeYerUqcyaNYvFixezatUqBgwYwIgRI9hrr704/PDDN9t0003vKdxuxowZ/OQnP8HdaWtrY/Dgwdu899577f369ePggw/mN7/5DQAzZ878yGD3oUOHsmTJErbddlujxPc1fvx4HnzwwdzTomUuvPBCzjrrLABuuOGGUjPP1CWX90888QTTpk1jzpw5vPrqq6xZs4aBAwey/fbbs++++3LwwQeP7dOnz6LC7efOnct1113HtGnTWLJkCeuvvz7jxo07mnC3hqo+g/y2tLa2cscdd3DPPffw7LPPsmLFCvr06cPWW2/N7rvvzoQJE/ptvfXWU4ApperZY489vlBmX2XH3MWxmv9iZjcB/0qYs3VTwsTmdxPmOM0
"text/plain": [
"<Figure size 480x300 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"density = True\n",
"cumulative = True\n",
"histtype = 'step'\n",
"lw = 2\n",
"if perform_zscore:\n",
" bins = {\n",
" 'stim_energy': np.arange(0, .7, .01),\n",
" 'stim_half_width': np.arange(0, 10, .5),\n",
" 'stim_p_max': np.arange(0, 4, .01),\n",
" 'stim_strength': np.arange(0, 160, 1)\n",
" }\n",
"else:\n",
" bins = {\n",
" 'stim_energy': np.arange(0, .008, .0001),\n",
" 'stim_half_width': np.arange(0, 0.5, .001),\n",
" 'stim_p_max': np.arange(0, .06, .0001),\n",
" 'stim_strength': np.arange(0, 160, 1)\n",
" }\n",
"xlabel = {\n",
" 'stim_energy': 'Energy (dB)',\n",
" 'stim_half_width': '(Hz)',\n",
" 'stim_p_max': 'Peak PSD (dB/Hz)',\n",
" 'stim_strength': 'Ratio',\n",
"}\n",
"# key = 'theta_energy'\n",
"# key = 'theta_peak'\n",
"for key in bins:\n",
" results[key] = list()\n",
" fig = plt.figure(figsize=(3.2,2))\n",
" plt.suptitle(key)\n",
" legend_lines = []\n",
" for color, query, label in zip(colors[1::2], queries[1::2], labels[1::2]):\n",
" values = lfp_results_hemisphere.query(query).loc[:,['entity_date_side', key]]\n",
" results[key].append(values.rename({key: label}, axis=1))\n",
" values[key].hist(\n",
" bins=bins[key], density=density, cumulative=cumulative, lw=lw, \n",
" histtype=histtype, color=color)\n",
" legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, lw=lw, label=label))\n",
" \n",
" plt.legend(\n",
" handles=legend_lines,\n",
" bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)\n",
" plt.tight_layout()\n",
" plt.grid(False)\n",
" plt.xlim(right=bins[key].max() - bins[key].max()*0.025)\n",
" despine()\n",
" plt.xlabel(xlabel[key])\n",
" figname = f'lfp-psd-histogram-{key}'\n",
" fig.savefig(\n",
" output_path / 'figures' / f'{figname}.png', \n",
" bbox_inches='tight', transparent=True)\n",
" fig.savefig(\n",
" output_path / 'figures' / f'{figname}.svg', \n",
" bbox_inches='tight', transparent=True)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"from functools import reduce"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"for key, val in results.items():\n",
" df = reduce(lambda left,right: pd.merge(left, right, on='entity_date_side', how='outer'), val)\n",
" results[key] = df.drop('entity_date_side',axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"def summarize(data):\n",
" return \"{:.1e} ± {:.1e} ({})\".format(data.mean(), data.sem(), sum(~np.isnan(data)))\n",
"\n",
"\n",
"def MWU(df, keys):\n",
" '''\n",
" Mann Whitney U\n",
" '''\n",
" Uvalue, pvalue = scipy.stats.mannwhitneyu(\n",
" df[keys[0]].dropna(), \n",
" df[keys[1]].dropna(),\n",
" alternative='two-sided')\n",
"\n",
" return \"{:.2f}, {:.3f}\".format(Uvalue, pvalue)\n",
"\n",
"\n",
"def PRS(df, keys):\n",
" '''\n",
" Permutation ReSampling\n",
" '''\n",
" pvalue, observed_diff, diffs = permutation_resampling(\n",
" df[keys[0]].dropna(), \n",
" df[keys[1]].dropna(), statistic=np.median)\n",
"\n",
" return \"{:.2f}, {:.3f}\".format(observed_diff, pvalue)\n",
"\n",
"\n",
"def wilcoxon(df, keys):\n",
" dff = df.loc[:,keys].dropna()\n",
" if len(dff[keys].dropna()) == 0:\n",
" statistic, pvalue = np.nan, np.nan\n",
" else:\n",
" statistic, pvalue = scipy.stats.wilcoxon(\n",
" dff[keys[0]], \n",
" dff[keys[1]],\n",
" alternative='two-sided')\n",
"\n",
" return \"{:.2f}, {:.3f}, ({})\".format(statistic, pvalue, len(dff))\n",
"\n",
"\n",
"def paired_t(df, keys):\n",
" dff = df.loc[:,[keys[0], keys[1]]].dropna()\n",
" statistic, pvalue = scipy.stats.ttest_rel(\n",
" dff[keys[0]], \n",
" dff[keys[1]])\n",
"\n",
" return \"{:.2f}, {:.3f}\".format(statistic, pvalue)\n",
"\n",
" \n",
"def normality(df, key):\n",
" if len(df[key].dropna()) < 8:\n",
" statistic, pvalue = np.nan, np.nan\n",
" else:\n",
" statistic, pvalue = scipy.stats.normaltest(\n",
" df[key].dropna())\n",
"\n",
" return \"{:.2f}, {:.3f}\".format(statistic, pvalue)\n",
"\n",
"\n",
"def shapiro(df, key):\n",
" if len(df[key].dropna()) < 8:\n",
" statistic, pvalue = np.nan, np.nan\n",
" else:\n",
" statistic, pvalue = scipy.stats.shapiro(\n",
" df[key].dropna())\n",
"\n",
" return \"{:.2f}, {:.3f}\".format(statistic, pvalue)\n",
"\n",
"def rename(name):\n",
" return name.replace(\"_field\", \"-field\").replace(\"_\", \" \").capitalize()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'theta_energy': Baseline I 11 Hz Baseline II 30 Hz\n",
" 0 0.003500 NaN NaN NaN\n",
" 1 0.001237 NaN NaN NaN\n",
" 2 0.003269 NaN 0.004628 NaN\n",
" 3 0.001465 NaN 0.002155 NaN\n",
" 4 0.002685 NaN NaN NaN\n",
" 5 0.000772 NaN NaN NaN\n",
" 6 0.002735 NaN 0.002988 NaN\n",
" 7 0.000872 NaN 0.001177 NaN\n",
" 8 0.004156 NaN 0.003132 NaN\n",
" 9 0.001317 NaN 0.001268 NaN\n",
" 10 0.000423 NaN 0.003200 NaN\n",
" 11 0.000957 NaN 0.001153 NaN\n",
" 12 0.002135 NaN 0.001651 NaN\n",
" 13 0.001649 NaN 0.000871 NaN\n",
" 14 0.009256 NaN 0.006838 NaN\n",
" 15 0.003803 NaN 0.003271 NaN\n",
" 16 0.004999 NaN 0.003724 NaN\n",
" 17 0.001193 NaN 0.002239 NaN\n",
" 18 0.004548 0.001780 NaN NaN\n",
" 19 0.004548 0.001093 NaN NaN\n",
" 20 0.002307 0.001479 NaN NaN\n",
" 21 0.002307 0.001122 NaN NaN\n",
" 22 0.002577 NaN 0.002665 NaN\n",
" 23 0.001591 NaN 0.001562 NaN\n",
" 24 0.001850 NaN 0.002546 NaN\n",
" 25 0.001388 NaN 0.001375 NaN\n",
" 26 0.002535 NaN 0.003099 NaN\n",
" 27 0.001129 NaN 0.001638 NaN\n",
" 28 0.001581 NaN 0.001302 NaN\n",
" 29 0.001346 NaN 0.000808 NaN\n",
" 30 0.002590 NaN 0.001583 NaN\n",
" 31 0.001029 NaN 0.000921 NaN\n",
" 32 0.001280 NaN NaN NaN\n",
" 33 0.000644 NaN NaN NaN\n",
" 34 0.002180 NaN NaN NaN\n",
" 35 0.002834 NaN NaN NaN\n",
" 36 0.001891 0.001814 NaN NaN\n",
" 37 0.001891 0.001557 NaN NaN\n",
" 38 0.002443 0.003549 NaN NaN\n",
" 39 0.002443 0.002262 NaN NaN\n",
" 40 0.002762 NaN 0.001113 NaN\n",
" 41 0.002649 NaN 0.003236 NaN\n",
" 42 0.000662 NaN NaN NaN\n",
" 43 0.003067 NaN NaN NaN\n",
" 44 0.001097 NaN 0.000649 NaN\n",
" 45 0.002022 NaN 0.001838 NaN\n",
" 46 0.001199 NaN 0.001475 NaN\n",
" 47 0.004233 NaN 0.001936 NaN\n",
" 48 0.002117 NaN 0.004077 NaN\n",
" 49 0.002906 NaN 0.003994 NaN,\n",
" 'theta_peak': Baseline I 11 Hz Baseline II 30 Hz\n",
" 0 0.003070 NaN NaN NaN\n",
" 1 0.001036 NaN NaN NaN\n",
" 2 0.003435 NaN 0.003384 NaN\n",
" 3 0.001416 NaN 0.001389 NaN\n",
" 4 0.002043 NaN NaN NaN\n",
" 5 0.000472 NaN NaN NaN\n",
" 6 0.002429 NaN 0.003210 NaN\n",
" 7 0.000625 NaN 0.000961 NaN\n",
" 8 0.002038 NaN 0.002704 NaN\n",
" 9 0.000872 NaN 0.001010 NaN\n",
" 10 0.000342 NaN 0.002434 NaN\n",
" 11 0.000640 NaN 0.000825 NaN\n",
" 12 0.001660 NaN 0.001235 NaN\n",
" 13 0.001190 NaN 0.000620 NaN\n",
" 14 0.007286 NaN 0.006490 NaN\n",
" 15 0.002671 NaN 0.002261 NaN\n",
" 16 0.004189 NaN 0.003590 NaN\n",
" 17 0.000846 NaN 0.001383 NaN\n",
" 18 0.003873 0.001374 NaN NaN\n",
" 19 0.003873 0.000699 NaN NaN\n",
" 20 0.001684 0.001046 NaN NaN\n",
" 21 0.001684 0.000645 NaN NaN\n",
" 22 0.002021 NaN 0.001967 NaN\n",
" 23 0.001159 NaN 0.001171 NaN\n",
" 24 0.001540 NaN 0.002067 NaN\n",
" 25 0.001135 NaN 0.001147 NaN\n",
" 26 0.002134 NaN 0.002335 NaN\n",
" 27 0.000959 NaN 0.001139 NaN\n",
" 28 0.001189 NaN 0.000792 NaN\n",
" 29 0.000657 NaN 0.000317 NaN\n",
" 30 0.001134 NaN 0.001158 NaN\n",
" 31 0.000427 NaN 0.000362 NaN\n",
" 32 0.000670 NaN NaN NaN\n",
" 33 0.000301 NaN NaN NaN\n",
" 34 0.001377 NaN NaN NaN\n",
" 35 0.001044 NaN NaN NaN\n",
" 36 0.001698 0.001544 NaN NaN\n",
" 37 0.001698 0.000898 NaN NaN\n",
" 38 0.002605 0.002939 NaN NaN\n",
" 39 0.002605 0.001677 NaN NaN\n",
" 40 0.002632 NaN 0.000722 NaN\n",
" 41 0.002425 NaN 0.003022 NaN\n",
" 42 0.000412 NaN NaN NaN\n",
" 43 0.003115 NaN NaN NaN\n",
" 44 0.000959 NaN 0.000390 NaN\n",
" 45 0.002131 NaN 0.001653 NaN\n",
" 46 0.000515 NaN 0.000891 NaN\n",
" 47 0.002570 NaN 0.001760 NaN\n",
" 48 0.001693 NaN 0.002711 NaN\n",
" 49 0.003029 NaN 0.002933 NaN,\n",
" 'theta_freq': Baseline I 11 Hz Baseline II 30 Hz\n",
" 0 7.4 NaN NaN NaN\n",
" 1 7.4 NaN NaN NaN\n",
" 2 7.8 NaN 7.8 NaN\n",
" 3 7.8 NaN 7.8 NaN\n",
" 4 7.4 NaN NaN NaN\n",
" 5 7.4 NaN NaN NaN\n",
" 6 7.5 NaN 7.8 NaN\n",
" 7 7.5 NaN 7.8 NaN\n",
" 8 7.6 NaN 7.9 NaN\n",
" 9 7.6 NaN 7.9 NaN\n",
" 10 7.4 NaN 7.7 NaN\n",
" 11 7.4 NaN 7.7 NaN\n",
" 12 7.6 NaN 7.6 NaN\n",
" 13 7.5 NaN 7.6 NaN\n",
" 14 8.3 NaN 8.5 NaN\n",
" 15 8.2 NaN 8.5 NaN\n",
" 16 7.9 NaN 8.1 NaN\n",
" 17 7.7 NaN 8.2 NaN\n",
" 18 8.0 8.2 NaN NaN\n",
" 19 8.0 8.2 NaN NaN\n",
" 20 8.0 8.2 NaN NaN\n",
" 21 8.0 8.2 NaN NaN\n",
" 22 8.1 NaN 8.2 NaN\n",
" 23 8.1 NaN 8.2 NaN\n",
" 24 8.0 NaN 8.2 NaN\n",
" 25 8.1 NaN 8.2 NaN\n",
" 26 8.1 NaN 8.2 NaN\n",
" 27 8.1 NaN 8.2 NaN\n",
" 28 8.0 NaN 8.2 NaN\n",
" 29 6.3 NaN 8.3 NaN\n",
" 30 8.2 NaN 8.6 NaN\n",
" 31 6.2 NaN 8.3 NaN\n",
" 32 8.3 NaN NaN NaN\n",
" 33 8.4 NaN NaN NaN\n",
" 34 8.2 NaN NaN NaN\n",
" 35 7.6 NaN NaN NaN\n",
" 36 7.9 7.9 NaN NaN\n",
" 37 7.9 7.9 NaN NaN\n",
" 38 7.9 8.0 NaN NaN\n",
" 39 7.9 7.9 NaN NaN\n",
" 40 7.6 NaN 8.1 NaN\n",
" 41 7.6 NaN 8.1 NaN\n",
" 42 7.4 NaN NaN NaN\n",
" 43 7.6 NaN NaN NaN\n",
" 44 7.7 NaN 7.8 NaN\n",
" 45 7.7 NaN 7.8 NaN\n",
" 46 7.8 NaN 8.0 NaN\n",
" 47 8.0 NaN 8.0 NaN\n",
" 48 7.7 NaN 8.3 NaN\n",
" 49 7.7 NaN 8.0 NaN,\n",
" 'theta_half_width': Baseline I 11 Hz Baseline II 30 Hz\n",
" 0 0.823669 NaN NaN NaN\n",
" 1 0.845791 NaN NaN NaN\n",
" 2 0.577806 NaN 0.959215 NaN\n",
" 3 0.653362 NaN 1.056302 NaN\n",
" 4 1.059779 NaN NaN NaN\n",
" 5 1.147158 NaN NaN NaN\n",
" 6 0.787100 NaN 0.637116 NaN\n",
" 7 0.964045 NaN 0.633677 NaN\n",
" 8 1.285535 NaN 0.760866 NaN\n",
" 9 0.913884 NaN 0.766207 NaN\n",
" 10 0.776057 NaN 0.996439 NaN\n",
" 11 0.928426 NaN 0.941342 NaN\n",
" 12 0.966327 NaN 1.072370 NaN\n",
" 13 0.984655 NaN 1.106900 NaN\n",
" 14 0.960731 NaN 0.738833 NaN\n",
" 15 1.044175 NaN 0.879523 NaN\n",
" 16 0.911633 NaN 0.641326 NaN\n",
" 17 0.953731 NaN 0.744352 NaN\n",
" 18 0.749468 0.930533 NaN NaN\n",
" 19 0.749468 0.935671 NaN NaN\n",
" 20 0.930580 0.969152 NaN NaN\n",
" 21 0.930580 1.049631 NaN NaN\n",
" 22 1.009373 NaN 1.085564 NaN\n",
" 23 1.065057 NaN 0.944394 NaN\n",
" 24 0.831508 NaN 0.910906 NaN\n",
" 25 0.825350 NaN 0.831260 NaN\n",
" 26 0.844569 NaN 0.885471 NaN\n",
" 27 0.720270 NaN 0.915350 NaN\n",
" 28 0.873920 NaN 1.122807 NaN\n",
" 29 NaN NaN 1.657428 NaN\n",
" 30 1.224911 NaN 1.120201 NaN\n",
" 31 6.316287 NaN NaN NaN\n",
" 32 1.015663 NaN NaN NaN\n",
" 33 1.318961 NaN NaN NaN\n",
" 34 1.124338 NaN NaN NaN\n",
" 35 8.080786 NaN NaN NaN\n",
" 36 0.755843 0.620998 NaN NaN\n",
" 37 0.755843 1.215511 NaN NaN\n",
" 38 0.650319 0.620433 NaN NaN\n",
" 39 0.650319 0.745154 NaN NaN\n",
" 40 0.677682 NaN 1.000266 NaN\n",
" 41 0.679729 NaN 0.541226 NaN\n",
" 42 0.945508 NaN NaN NaN\n",
" 43 0.604237 NaN NaN NaN\n",
" 44 0.583246 NaN 0.870935 NaN\n",
" 45 0.560840 NaN 0.754164 NaN\n",
" 46 0.843262 NaN 1.160272 NaN\n",
" 47 0.969953 NaN 0.680846 NaN\n",
" 48 0.823915 NaN 0.979878 NaN\n",
" 49 0.688965 NaN 0.850279 NaN,\n",
" 'stim_energy': 11 Hz 30 Hz\n",
" 0 0.002671 NaN\n",
" 1 0.000532 NaN\n",
" 2 0.000168 NaN\n",
" 3 0.000033 NaN\n",
" 4 0.000513 NaN\n",
" 5 0.002245 NaN\n",
" 6 0.000340 NaN\n",
" 7 0.000240 NaN,\n",
" 'stim_half_width': 11 Hz 30 Hz\n",
" 0 10.318750 NaN\n",
" 1 0.149163 NaN\n",
" 2 0.156917 NaN\n",
" 3 0.271082 NaN\n",
" 4 0.150665 NaN\n",
" 5 6.113351 NaN\n",
" 6 0.655030 NaN\n",
" 7 0.156293 NaN,\n",
" 'stim_p_max': 11 Hz 30 Hz\n",
" 0 0.000105 NaN\n",
" 1 0.004951 NaN\n",
" 2 0.001566 NaN\n",
" 3 0.000206 NaN\n",
" 4 0.004785 NaN\n",
" 5 0.000374 NaN\n",
" 6 0.000814 NaN\n",
" 7 0.002239 NaN,\n",
" 'stim_strength': 11 Hz 30 Hz\n",
" 0 0.059252 NaN\n",
" 1 4.530157 NaN\n",
" 2 1.059344 NaN\n",
" 3 0.183414 NaN\n",
" 4 2.637518 NaN\n",
" 5 0.239875 NaN\n",
" 6 0.229359 NaN\n",
" 7 0.989594 NaN}"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/stats/morestats.py:2863: UserWarning: Sample size too small for normal approximation.\n",
" warnings.warn(\"Sample size too small for normal approximation.\")\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/stats/stats.py:1450: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=8\n",
" \"anyway, n=%i\" % int(n))\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Theta energy</th>\n",
" <th>Theta peak</th>\n",
" <th>Theta freq</th>\n",
" <th>Theta half width</th>\n",
" <th>Stim energy</th>\n",
" <th>Stim half width</th>\n",
" <th>Stim p max</th>\n",
" <th>Stim strength</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>11 Hz</th>\n",
" <td>1.8e-03 ± 2.8e-04 (8)</td>\n",
" <td>1.4e-03 ± 2.6e-04 (8)</td>\n",
" <td>8.1e+00 ± 5.3e-02 (8)</td>\n",
" <td>8.9e-01 ± 7.4e-02 (8)</td>\n",
" <td>8.4e-04 ± 3.6e-04 (8)</td>\n",
" <td>2.2e+00 ± 1.4e+00 (8)</td>\n",
" <td>1.9e-03 ± 7.0e-04 (8)</td>\n",
" <td>1.2e+00 ± 5.6e-01 (8)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>30 Hz</th>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" <td>nan ± nan (0)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Baseline I</th>\n",
" <td>2.3e-03 ± 2.1e-04 (50)</td>\n",
" <td>1.8e-03 ± 1.8e-04 (50)</td>\n",
" <td>7.8e+00 ± 5.9e-02 (50)</td>\n",
" <td>1.1e+00 ± 1.8e-01 (49)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Baseline II</th>\n",
" <td>2.3e-03 ± 2.4e-04 (32)</td>\n",
" <td>1.8e-03 ± 2.3e-04 (32)</td>\n",
" <td>8.1e+00 ± 4.7e-02 (32)</td>\n",
" <td>9.1e-01 ± 3.9e-02 (31)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Normality 11 Hz</th>\n",
" <td>8.13, 0.017</td>\n",
" <td>6.61, 0.037</td>\n",
" <td>6.30, 0.043</td>\n",
" <td>0.23, 0.890</td>\n",
" <td>3.45, 0.178</td>\n",
" <td>7.42, 0.024</td>\n",
" <td>1.81, 0.405</td>\n",
" <td>6.43, 0.040</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Normality 30 Hz</th>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" <td>nan, nan</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Normality Baseline I</th>\n",
" <td>41.72, 0.000</td>\n",
" <td>31.09, 0.000</td>\n",
" <td>29.47, 0.000</td>\n",
" <td>81.74, 0.000</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Normality Baseline II</th>\n",
" <td>13.17, 0.001</td>\n",
" <td>20.78, 0.000</td>\n",
" <td>0.96, 0.618</td>\n",
" <td>13.33, 0.001</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon 11 Hz 30 Hz</th>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon 11 Hz Baseline II</th>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon Baseline I 11 Hz</th>\n",
" <td>5.00, 0.069, (8)</td>\n",
" <td>2.00, 0.025, (8)</td>\n",
" <td>0.00, 0.034, (8)</td>\n",
" <td>6.00, 0.093, (8)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon Baseline I 30 Hz</th>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon Baseline I Baseline II</th>\n",
" <td>264.00, 1.000, (32)</td>\n",
" <td>256.00, 0.881, (32)</td>\n",
" <td>0.00, 0.000, (32)</td>\n",
" <td>203.00, 0.544, (30)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Wilcoxon Baseline II 30 Hz</th>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>nan, nan, (0)</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Theta energy \\\n",
"11 Hz 1.8e-03 ± 2.8e-04 (8) \n",
"30 Hz nan ± nan (0) \n",
"Baseline I 2.3e-03 ± 2.1e-04 (50) \n",
"Baseline II 2.3e-03 ± 2.4e-04 (32) \n",
"Normality 11 Hz 8.13, 0.017 \n",
"Normality 30 Hz nan, nan \n",
"Normality Baseline I 41.72, 0.000 \n",
"Normality Baseline II 13.17, 0.001 \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II nan, nan, (0) \n",
"Wilcoxon Baseline I 11 Hz 5.00, 0.069, (8) \n",
"Wilcoxon Baseline I 30 Hz nan, nan, (0) \n",
"Wilcoxon Baseline I Baseline II 264.00, 1.000, (32) \n",
"Wilcoxon Baseline II 30 Hz nan, nan, (0) \n",
"\n",
" Theta peak \\\n",
"11 Hz 1.4e-03 ± 2.6e-04 (8) \n",
"30 Hz nan ± nan (0) \n",
"Baseline I 1.8e-03 ± 1.8e-04 (50) \n",
"Baseline II 1.8e-03 ± 2.3e-04 (32) \n",
"Normality 11 Hz 6.61, 0.037 \n",
"Normality 30 Hz nan, nan \n",
"Normality Baseline I 31.09, 0.000 \n",
"Normality Baseline II 20.78, 0.000 \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II nan, nan, (0) \n",
"Wilcoxon Baseline I 11 Hz 2.00, 0.025, (8) \n",
"Wilcoxon Baseline I 30 Hz nan, nan, (0) \n",
"Wilcoxon Baseline I Baseline II 256.00, 0.881, (32) \n",
"Wilcoxon Baseline II 30 Hz nan, nan, (0) \n",
"\n",
" Theta freq \\\n",
"11 Hz 8.1e+00 ± 5.3e-02 (8) \n",
"30 Hz nan ± nan (0) \n",
"Baseline I 7.8e+00 ± 5.9e-02 (50) \n",
"Baseline II 8.1e+00 ± 4.7e-02 (32) \n",
"Normality 11 Hz 6.30, 0.043 \n",
"Normality 30 Hz nan, nan \n",
"Normality Baseline I 29.47, 0.000 \n",
"Normality Baseline II 0.96, 0.618 \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II nan, nan, (0) \n",
"Wilcoxon Baseline I 11 Hz 0.00, 0.034, (8) \n",
"Wilcoxon Baseline I 30 Hz nan, nan, (0) \n",
"Wilcoxon Baseline I Baseline II 0.00, 0.000, (32) \n",
"Wilcoxon Baseline II 30 Hz nan, nan, (0) \n",
"\n",
" Theta half width \\\n",
"11 Hz 8.9e-01 ± 7.4e-02 (8) \n",
"30 Hz nan ± nan (0) \n",
"Baseline I 1.1e+00 ± 1.8e-01 (49) \n",
"Baseline II 9.1e-01 ± 3.9e-02 (31) \n",
"Normality 11 Hz 0.23, 0.890 \n",
"Normality 30 Hz nan, nan \n",
"Normality Baseline I 81.74, 0.000 \n",
"Normality Baseline II 13.33, 0.001 \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II nan, nan, (0) \n",
"Wilcoxon Baseline I 11 Hz 6.00, 0.093, (8) \n",
"Wilcoxon Baseline I 30 Hz nan, nan, (0) \n",
"Wilcoxon Baseline I Baseline II 203.00, 0.544, (30) \n",
"Wilcoxon Baseline II 30 Hz nan, nan, (0) \n",
"\n",
" Stim energy Stim half width \\\n",
"11 Hz 8.4e-04 ± 3.6e-04 (8) 2.2e+00 ± 1.4e+00 (8) \n",
"30 Hz nan ± nan (0) nan ± nan (0) \n",
"Baseline I NaN NaN \n",
"Baseline II NaN NaN \n",
"Normality 11 Hz 3.45, 0.178 7.42, 0.024 \n",
"Normality 30 Hz nan, nan nan, nan \n",
"Normality Baseline I NaN NaN \n",
"Normality Baseline II NaN NaN \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II NaN NaN \n",
"Wilcoxon Baseline I 11 Hz NaN NaN \n",
"Wilcoxon Baseline I 30 Hz NaN NaN \n",
"Wilcoxon Baseline I Baseline II NaN NaN \n",
"Wilcoxon Baseline II 30 Hz NaN NaN \n",
"\n",
" Stim p max Stim strength \n",
"11 Hz 1.9e-03 ± 7.0e-04 (8) 1.2e+00 ± 5.6e-01 (8) \n",
"30 Hz nan ± nan (0) nan ± nan (0) \n",
"Baseline I NaN NaN \n",
"Baseline II NaN NaN \n",
"Normality 11 Hz 1.81, 0.405 6.43, 0.040 \n",
"Normality 30 Hz nan, nan nan, nan \n",
"Normality Baseline I NaN NaN \n",
"Normality Baseline II NaN NaN \n",
"Wilcoxon 11 Hz 30 Hz nan, nan, (0) nan, nan, (0) \n",
"Wilcoxon 11 Hz Baseline II NaN NaN \n",
"Wilcoxon Baseline I 11 Hz NaN NaN \n",
"Wilcoxon Baseline I 30 Hz NaN NaN \n",
"Wilcoxon Baseline I Baseline II NaN NaN \n",
"Wilcoxon Baseline II 30 Hz NaN NaN "
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"\n",
"stat = pd.DataFrame()\n",
"\n",
"for key, df in results.items():\n",
" Key = rename(key)\n",
" stat[Key] = df.agg(summarize)\n",
" stat[Key] = df.agg(summarize)\n",
"\n",
" for i, c1 in enumerate(df.columns):\n",
" stat.loc[f'Normality {c1}', Key] = normality(df, c1)\n",
"# stat.loc[f'Shapiro {c1}', Key] = shapiro(df, c1)\n",
" for c2 in df.columns[i+1:]:\n",
"# stat.loc[f'MWU {c1} {c2}', Key] = MWU(df, [c1, c2])\n",
"# stat.loc[f'PRS {c1} {c2}', Key] = PRS(df, [c1, c2])\n",
" stat.loc[f'Wilcoxon {c1} {c2}', Key] = wilcoxon(df, [c1, c2])\n",
"# stat.loc[f'Paired T {c1} {c2}', Key] = paired_t(df, [c1, c2])\n",
"\n",
"stat.sort_index()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"stat.to_latex(output_path / \"statistics\" / \"statistics.tex\")\n",
"stat.to_csv(output_path / \"statistics\" / \"statistics.csv\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Plot PSD"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"psd = pd.read_feather(pathlib.Path(\"output\") / (\"stimulus-lfp-response\" + zscore_str) / 'data' / 'psd.feather')\n",
"freqs = pd.read_feather(pathlib.Path(\"output\") / (\"stimulus-lfp-response\" + zscore_str) / 'data' / 'freqs.feather')"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"from septum_mec.analysis.plotting import plot_bootstrap_timeseries"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"freq = freqs.T.iloc[0].values\n",
"\n",
"mask = (freq < 49)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3257: RuntimeWarning: Mean of empty slice.\n",
" out=out, **kwargs)\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
" ret = ret.dtype.type(ret / rcount)\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/_methods.py:154: RuntimeWarning: invalid value encountered in true_divide\n",
" ret, rcount, out=ret, casting='unsafe', subok=False)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAFFCAYAAADsEyV2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZxcVZnw8d+9tVf13umkk85KQm4IkSiIrCIhIhJGRARHUXwFAVl0hE9cMggOg4KADC6oIMwg66iAYQkOOspmAk4iSgghySVbd9Kd7vRee939/eNWV6rS3UnvS+V8+YTurrrLqerqU0+d+5znSI7jIAiCIAiCIAjFSB7vBgiCIAiCIAjCaBHBriAIgiAIglC0RLArCIIgCIIgFC0R7AqCIAiCIAhFSwS7giAIgiAIQtESwa4gCIIgCIJQtESwKwiCIAiCIBQtEewKgiAIgiAIRUsEu4IgCIIgCELREsGuIAiCIAiCULREsCsIgiAIgiAULRHsCoIgCIIgCEVLBLuCIAiCIAhC0RLB7gSlKMrziqI8P97tEARBGAuizxMEYbR4x7sBQr/mL1iwYDHgjHdDBEGYsKTxbsAIEn2eIAiHM6Q+T4zsCoIgCIIgCEVLBLuCIAiCIAhC0RLBriAIgiAIglC0RLArCIIgCIIgFC0R7AqCIAiCIAhFSwS7giAIgiAIQtESwa4gCIIgCIJQtESwKwiCIAiCIBQtEewKE56V6CC56Q/Y6fh4N0UQBEEQhElGrKAmTHjR1/4LW09hdu+j/IzLx7s5giAIgiBMImJkV5jwbD0NjoPZ1TTeTREEQRAEYZIRwa4wKTi2DY4z3s0QBEEQBGGSEcGuMDk49ni3QBAEQRCESUgEu8Lk4Ng4YmRXEARBEIRBEsGuMDmIkV1BEARBEIZABLvCpOCIYFcQBEEQhCEQwa4wOYhgVxAEQRCEIRDBrjA52CLYFQRBEARh8ESwK0wKIo1BEARBEIShEMGuMDmISgyCIAiCIAyBCHaFyUGM7AqCIAiCMAQi2BXGhNHZ6C77O1Qi2BUEQRAEYQhEsCuMOqOzkeirDxJ97T+HfAxHTFADoLGxEUVR+vy3aNEiTjrpJC644ALuueceOjs7x7u5Q7J69WoUReGMM84ouP3SSy9FURR+9KMfjVPLDm/VqlUoisKll1463k0RhElD9GuTt1+bLH2ed7wbIBQ/Y/8OrGQXOA52JoEcLBn8QcTIbi8LFy6kpOTAc2lZFtFolO3bt7N161aefPJJHnnkERRFGcdWCoIgDJzo14TRIIJdYfTJMo6RxkpLGG27Ccx63+CPIYLdXm666SZOOumkXrd3dXWxatUqXn31Vf7lX/6FF198EVme/Bdx7rzzTtLpNJWVlePdFEEQRono14TRMPlfKcKE52hJHFPHTscwu5qGdgwR7A5YZWUld9xxB36/n/r6etatWzfmbTDaG+h+5ZcYHXtH7JgzZsxg/vz5VFVVjdgxBUGYHCZCvzYaRL82NkSwK4w6W0viWAZ2Jo4ZbRnaQUSwOyiVlZUcffTRAGzfvn3Mzx9d+zBmVzOxtb8a83MLwnizbYddOzvo7EyNd1OKynj3a8LkJdIYhFFnZxI4pu7+04bY+Ys6u4NmmiYAkUikz/teeOEF/vCHP/Duu+/S3d2N1+tl6tSpnHTSSVx22WXMmzev135r167liSee4O233yYWi1FSUsLChQv5+Mc/zsUXX4zf789t6zgWOBKJRIJHHnmEP/3pTzQ0NOA4DrNmzeLss8/mS1/6EmVlZQN6PJdeeikbNmzg6quv5oYbbgDciS3Lly9nypQprFu3jqeffponn3ySHTt2AG7+32c+8xkuvPBCJEnqdcyRapsg5OvoSLJtSyser8ypp8+ltDQw3k0qGuPdr/UQ/drkIoJdYdTFkl1g6kiWia0fOth1HIe0aRTclnYcDNvBdBxShj6aTR1RIa+vz45oLOzZs4ft27cjyzIf/vCHC+7LZDJcddVVrF+/HoC6ujoWLlxIR0cH9fX11NfXs2bNGp544gkWL16c2+/RRx/ltttuA2Dq1KksWrSIrq4uNmzYwIYNG/jDH/7Aww8/jMfjcXewLRrak3zz/PNpamrC4/Ewa9YsgsEgO3bs4Oc//znPPvssDz74IPPnzx/W43Uch29/+9s899xzlJWVMW/ePPbu3cvGjRvZuHEju3fv5hvf+EbBPjt37uTKK68c9bYJR55ETCOV0tE1i3ffaeGkU2Zj6NZ4N2tE+PyeI7tfY+z6DtGvjZwjPthVFOUi4BrgBCAA7AVeAO5WVXXfIfY7DVgFnApEsvutBn6gqmr3aLd7smhORvlbo8rMeAe1gRCOnsJxnD47S8dx+NT/3M+brQ29DxSYAdEMPP7dMWj1yDhx6hxWr7h6zN4YLMsiFovx1ltvceedd2LbNldffTV1dXUF2z344IOsX7+eyspKHnjgAY477rjcfZs2beLaa6+lra2N+++/n5/+9KcAxGIx7r77bgDuuecezjvvvNw+69at47rrrsu9MfTcl87ofOuR12jqSLB8+XL+7d/+jWnTpgHQ1tbGTTfdxKuvvsq1117Lc889RzAYHPJj7+jo4IUXXuA73/kOn//85/F4PGiaxk033cTzzz/Pr371Ky6//PJcXlwqleKaa66hqalp1NsmHHkSSR1dt4jFMnR2prj3R+to3Bsd72aNiLnzKrnma6cesf3aWPYdol8bOUd0zq6iKP8JPAWcBaSBbcAM4Abg3WxA29d+nwH+AvxTdr93gVnAt4C3FUWZNfqtnxzq4x1kMnEyepqMnsGxLZxDLC4xPuMFk9MXv/jFgnqUixcv5uSTT+aaa66hvr6eK6+8kuuvv77Xfm+88QayLPPVr3614A0B4LjjjuNzn/scAO+9917u9t27d6NpGuXl5axYsaJgn9NPP52rrrqKc845B5/PB4AkSax5czeNHQmOPfZY7r333lynC1BTU8NPfvIT6urqqK+vZ/Xq1cN+Pi655BK++MUv5kZgAoEAN954I5IkYZommzZtym371FNP0dDQMGZtE44syYSOoVs4QCZtYNsiDWugJnK/BmPfd4h+bWQcsSO7iqJ8GfgyYAJfUlX1ieztZcADwD8DqxVFmaeqaipvPwV4DPeDwteAn6uq6iiKMg03cP4w8N/Zr0e89mQU09AwTZ2Ux0ulZeDoKQiEe20rSRKrV1zdK42hY83t6M3b8E9fRPUnbhyrpg/baKcxHFyP0rZtkskkDQ0N6LrOww8/TDqd5sYbbyy4BPfrX/8awzD6bVsoFALcy4I9Zs6cidfrJRqNsmrVKi677DIWLVqUu/+6664rPIjs4S9b3QsjK1asKDh/j2AwyDnnnMNDDz3EK6+8wiWXXDL4JyHPsmXLet1WWVlJVVUVHR0dxGKx3O1//vOfh9w2x9QBCcnr67WfIABomolp2gQCXkzT5rzzj2HOnOIoLTXaaQwTul9jeH3HUIxVv1bsjthgF/hm9usPewJdAFVVY4qiXAZ8DJgKfAp4Im+/fwX8wG9UVf1Z3n77FUW5ANgFnK4oykdVVf3zaD+Iia4j3o4v1Y0vHSUdiLglyLQkntIpfW4vSRJhX+FkgJQk4XEs/H3cdyTrrx6lrus888wzfO973+Pxxx/HsixuueWWgm18Ph/RaJSNGzdSX1/P3r17qa+vZ+vWrbS3twPum0yP6upqrrjiCu6//36effZZnn32WWpqajj55JM5/fTTOeOMMwpK50iyl92tccAdbXjppZf6fAw959q1a9ewngugYBQjX8+lOss6kDPZM7oz2LY5jkPXH36EY2pUnf8dJLn3G4pwZLNtB12zME2bsvIApmGTyZj4A0fy2+3ATeR+DYbedwzVWPRrR4Ij8q9PURQ/8DywhMJAFgBVVdOKomwHPgTMydsviDviC/BffezXqSjKU8AVwOeAIzrYNW2LRDrGSY1vUZ7uJubxYegZ7Ex80McSywUPnN/v55//+Z9pa2vj3nvv5cknn+QrX/kK06dPB9yZurfddhtr1qzBMA6Movt8Po499li
"text/plain": [
"<Figure size 750x300 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(5,2))\n",
"axs = axs.repeat(2)\n",
"for i, (ax, query) in enumerate(zip(axs.ravel(), queries)):\n",
" selection = [\n",
" f'{r.action}_{r.channel_group}' \n",
" for i, r in lfp_results_hemisphere.query(query).iterrows()]\n",
" values = psd.loc[mask, selection].to_numpy()\n",
" values = 10 * np.log10(values)\n",
" plot_bootstrap_timeseries(freq[mask], values, ax=ax, lw=1, label=labels[i], color=colors[i])\n",
"# ax.set_title(titles[i])\n",
" ax.set_xlabel('Frequency Hz')\n",
" ax.legend(frameon=False)\n",
"axs[0].set_ylabel('PSD (dB/Hz)')\n",
"# axs[0].set_ylim(-31, 1)\n",
"despine()\n",
"\n",
"figname = 'lfp-psd'\n",
"fig.savefig(\n",
" output_path / 'figures' / f'{figname}.png', \n",
" bbox_inches='tight', transparent=True)\n",
"fig.savefig(\n",
" output_path / 'figures' / f'{figname}.svg', \n",
" bbox_inches='tight', transparent=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Store results in Expipe action"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"action = project.require_action(\"stimulus-lfp-response-mec\" + zscore_str)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/statistics/statistics.tex',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/statistics/statistics.csv',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_energy.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_strength.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_peak.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_p_max.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_freq.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_energy.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_freq.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_half_width.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_half_width.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_half_width.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_energy.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_peak.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_p_max.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-theta_half_width.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/.ipynb_checkpoints/lfp-psd-histogram-theta_energy-checkpoint.png',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_energy.svg',\n",
" '/media/storage/expipe/septum-mec/actions/stimulus-lfp-response-mec-no-zscore/data/figures/lfp-psd-histogram-stim_strength.svg']"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"copy_tree(output_path, str(action.data_path()))"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"septum_mec.analysis.registration.store_notebook(action, \"20_stimulus-lfp-response-mec.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": 4
}