notebooks/09_IRASA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aperiodic and oscillatory components of the power spectrum\n",
    "\n",
    "This notebook demonstrates how to use YASA to separate the aperiodic (= fractal = $1/f$) components of the EEG power spectra using the IRASA method.\n",
    "\n",
    "Please make sure to install the latest version of YASA first by typing the following line in your terminal or command prompt:\n",
    "\n",
    "`pip install --upgrade yasa`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "import yasa\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(style='white', font_scale=1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data loading\n",
    "\n",
    "Let's load 6 min of eyes open EEG data, sampled at 200 Hz. We use the [MNE library](https://mne.tools/stable/index.html), which includes various functions and utilities for reading EEG data and electrode locations. For a description of all the file formats supported by MNE, please refer to https://mne.tools/stable/auto_tutorials/io/plot_20_reading_eeg_data.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chan = ['F4-A1', 'CZ-A2']\n",
      "Sampling frequency = 200.0 Hz\n",
      "Data shape = (2, 72000)\n",
      "Duration = 360.0 seconds\n"
     ]
    }
   ],
   "source": [
    "# Load data as a MNE Raw file\n",
    "raw = mne.io.read_raw_fif('data_resting_EO_200Hz_raw.fif', preload=True, verbose=0)\n",
    "\n",
    "# Extract data, sf, and chan\n",
    "data = raw.get_data(units=\"uV\")\n",
    "sf = raw.info['sfreq']\n",
    "chan = raw.ch_names\n",
    "\n",
    "# Let's have a look at the data\n",
    "print('Chan =', chan)\n",
    "print('Sampling frequency =', sf, 'Hz')\n",
    "print('Data shape =', data.shape)\n",
    "print('Duration =', data.shape[1] / sf, 'seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**************\n",
    "\n",
    "## Calculate the original power spectrum\n",
    "\n",
    "For a primer on EEG spectral power please refer to https://raphaelvallat.com/bandpower.html or [08_bandpower.ipynb](08_bandpower.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(401,) (2, 401)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEgCAYAAAD2c3e8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAG0lEQVR4nO3deViU9f7/8ecw7CAIiCCKsmgorkgBJmqKCyfNytL8uov+zDQ9uZbHUrNFSxMz9ZtI5XEvFDPTyo64a6S5looI5oILKotswwAzvz/8NqcJTZaBG8b347q4Ludzb++bW3hxb5+PSq/X6xFCCCEUYqF0AUIIIR5tEkRCCCEUJUEkhBBCURJEQgghFCVBJIQQQlESREIIIRQlQSRENcnNzSU6OppevXrRpk0bunXrxgcffEBOTg4AQ4cOJSAg4L5fb7zxxkPXf+TIEQICApg1a9Z9p+/evZvnn3+eoKAgIiMjiYuLM+n+CVFRKnmPSIiql5OTw8CBA3F0dOTVV1/Fx8eHixcv8uGHH2Jvb8+aNWsoKCigqKjIaLmdO3cyd+5cPv74YyIjI/92G2+99RY///wzd+7c4cCBA9ja2hqmnT59moEDBzJ16lQiIiI4duwYb775Jh999BG9evWqkn0WoqzkjEiIarBw4UL0ej2rVq2iU6dOeHt707lzZ2JiYjhz5gxbtmyhbt26uLu7G75ycnJYuHAhL7744kNDSKvV8sMPPzB27FgKCwv54YcfjKZ//fXXhIWFMXLkSBo3bsxzzz3Hc889x9atW6tyt4UoE0ulCxDC3Gm1WrZt28bUqVOxs7Mzmubl5cW///1v/P39Sy0zZcoUPDw8mDlz5kO3sWfPHnJycujSpQthYWFs3ryZZ5991jD9pZdeQqfTGS2jUqkMlwWFUJKcEQlRxa5cuUJeXh5t2rS57/Tg4GDq1q1r1BYdHU1ycjKLFi3C3t7+odv45ptvaN++Pa6urvTs2ZOff/6ZK1euGKY/9thjNG/e3PA5PT2d7du306FDh4rtlBAmJEEkRBXLzs4GoE6dOmWa/9ChQ3zxxRdMmTKFwMDAh85/9+5d9u7dS48ePQCIiIhArVazZcuW+86fm5vL+PHjcXd3Z/jw4WXcCyGqjgSREFXMxcUF+G8g/Z2MjAymT59Ox44dGTFihNG0b775hqCgIMPXp59+CsB3332HVqulZ8+eALi6uvLEE0/w9ddfl7ocl5mZyciRI7l+/ToxMTE4ODiYYA+FqBy5RyREFWvcuDF169bl9OnT9708N2/ePBo0aMCIESOYOXMmOp2ODz74AJVKZTRft27daNu2reGzs7MzcC+gALp3726YptPp0Ov1HD58mI4dOwJw8+ZNRo4cSUFBAevWraNJkyYm31chKkKCSIgqplar6dOnD2vWrOGFF14weqz60qVLbNiwgddff51169axe/duVqxYQb169Uqtx9HREUdHR6O2tLQ0fvnlFyZMmGA4IwIoLi5myJAhbN68mY4dO5Kbm0tUVBRarZa1a9fSsGHDqtthIcpJgkiIajB+/Hj27dvHiBEjmDBhAo0bN+bs2bN8+OGHtGzZkuDgYF566SX69etHYGAgt27dMlreysqq1AMNANu2bcPGxoZhw4bh5ORkNO35558nLi6O7Oxsli9fzpUrV1i1ahXW1taG9VtaWhouHQqhFHmhVYhqcvv2bZYvX86ePXu4ffs2Hh4e9OrVi7Fjx7Jz505mzJjxwGVDQkJYs2ZNqfbevXvTpk0b5s2bV2paamoqTz/9NG+99RYrVqzg5s2bpeZp2bIl8fHxldsxISpJgkgIIYSizPLSXHZ2NlOmTKGwsBBra2sWLlwolx+EEKKGMsvHt9euXUuvXr1Ys2YNkZGRrF27VumShBBCPIBZnhENGzYMa2trAIqKirCxsVG4IiGEEA9ilkH0xxvs586dY8OGDXJGJIQQNZhZXpoDOHbsGNOmTWPJkiWGF//Kori4mKtXr1JcXFyF1QkhhPiDWQZRcnIyc+bMYcWKFfj6+pZr2Rs3bhAREcGNGzeqqDohhBB/ZpaX5j755BM0Gg2vv/46cK9rlJEjRypclRBCiPsxyyBasmSJ0iUIIYQoo1pzae7UqVOEhIQYtSUlJTFw4EDatWtHZGQke/fuVag6IYQQFVUrgmjHjh1ERUVRVFRkaNNqtbzyyiv06tWLI0eOMG3aNF577TXS0tIUrFQIIUR51fggio6OJjY2lnHjxhm1JyYmotFoGDFiBFZWVkRERBASEsK2bdsUqlQIIURF1PggGjRoEPHx8bRs2dKoPSUlBX9/f6MxW/z8/EhOTq7uEoUQQlRCjQ8iDw+P+7bn5+cbjesCYGtrS0FBQXWUJYQQwkRqfBA9iL29PRqNxqhNo9Fgb2+vUEVCCCEqotYGkb+/PxcvXjRqS01NpWnTpgpVJIQQoiJqbRCFhoaiVquJiYlBq9WSkJBAYmIivXv3Vro0IYQQ5VBrg8ja2pqVK1eyZ88ewsLCWLBgAdHR0Xh7eytdmhBCiHKQEVr/4urVq0RERLBr1y4aNWqkdDlCCGH2au0ZkRBCCPMgQSSEEEJREkRCCCEUZZa9b4tHQ2JiIt9++y0uLi64u7vTpk0b2rZtq3RZQohykiAStVJxcTHjx4/n+vXrhja1Wk1CQoK8SyZELSOX5kSttHv3bq5fv45jXXcatPwHLp5NKSkpkbGohKiFJIhErbRu3ToA6vqE06D1M3g9PhwLtZotW7aU6nFDCFGzSRCJWuf69evs2rULtdoSF58wAGwc3fDwD0On07F06VKFKxRClIcEkah1Nm7ciE6nw92nHVa2dQztLs16oLKwYNOmTVy5ckXBCoUQ5SFBJGqVkpISNmzYAICjdwejabZ16uPpH0pxcTHLli1TojwhRAVIEIlaZd++faSlpeHk4kEdj4BS012a9UClUrFx40ZSUlIUqFAIUV4SRKJWWbNmDQBOjcNQqUr/97V18sSzWUeKiop4++23q7s8IUQFSBCJWuPy5cvs3LkTtaUVLr4dHjhfvcA+2Njas2vXLhISEqqxQiFERUgQiVrjiy++QK/XU9/3CaxsnR44n5WtE/UD/wHA7Nmz0Wq11VWiEKICJIge4M6dO0qXIP4kLy+PjRs3AuDs2/mh87v5d8G5nhepqal8/vnnVV2eEKISJIge4Msvv1S6BPEnmzZt4u7du9Rr+Bj2ro0fOr+F2pJ6gc8DsHz5coqLi6u6RCFEBUkQPcCmTZvIzMxUugwB6PV6w1lNnSadyrycU4NAnN0acOfOHY4cOVJV5QkhKkmC6AHy8/OJjY1VugzBvUe2L1y4gKOzG3UbtSvzciqVCkfPNgDs2LGjiqoTQlSWBNHf+Oyzz8jOzla6jEeaTqfjww8/BO71K6eyUJdreUeve8NC7NixA51OZ/L6hBCVJ0H0AK4NHiMnJ4fPPvtM6VIeaXFxcZw4cQIHJ1fcmnYp9/L2ro1xdK7HjRs3OHnyZBVUKISoLLMPojVr1hATE1Pu5Zx8792LiI2NlXtFCrl79y7vv/8+AO6BfVFb2ZZ7HSqVBU5e9y7PfffddyatTwhhGmYdRHPmzGH16tUVWtbRzYd6jVqQnZ3NwIEDuX37tomrEw+zePFibt++jZtXM1yaPFHh9Tj8332i7du3o9frTVWeEMJEzDqIwsPDeeWVVyq8vGfQYJxcPfn11195/vnnSUtLM2F14u9cuHCBzz77DJVKhXurF1CpVBVel6O7P3YOTvz+++8kJSWZsEohhCmYdRB17969UstbO7jSpPNruNRvQmpqKs8995yEUTX5+OOPKS4uxvOx8DK9N/R3VBZqXBrde2hBLs8JUfOYdRCZgpWtE97hE3Dzasa1a9cYPnw4OTk5Spdl1q5cucLWrVuxsFDj2qynSdZp59EauDeya35+vknWKYQwDQmiMrC0tqdh2Ms4u3lx9uxZxo0bJ2/qV6GYmBhKSkrw8A/BxtHNJOt08gzE1dOX69ev8/HHH5tknUII05AgKqN7YTQGOwcnEhISZIiBKpKRkcH69esBqOvfzWTrVVlY4N76RQBWrFjBhQsXTLZuIUTlmH0Q9evXjzFjxphkXbZ16uMdOhpLSys+//xzvv32W5OsV/zX559/jkajwcOnLXZ1G5p03Q5uvjQI6ERRUREzZ86UJ+iEqCFqVRCdOnWKkJAQo7akpCQGDhxIu3btiIyMZO/evVVag2P9pni16wfA66+/zo0bN6p0e4+S/Px8vvjiCwCc/SKqZBvugc9gZ1+HAwcO8M0331TJNoQQ5VNrgmjHjh1ERUVRVFRkaNNqtbzyyiv06tWLI0eOMG3aNF577bUqf7LNzb8z7o3bkJWVxeTJk+UvaxOJjY0lKyuLel7NcHD3r5JtWNo4Ur/lMwDMnDlT/pAQogaoFUEUHR1NbGws48aNM2pPTExEo9EwYsQIrKysiIiIICQkhG3btlVpPSqVCs+g/8HOwYm9e/eyatWqKt3eo+D27dssW7YMAJeApyv13tDDuPo9iXvjVmRmZjJ58mTpg04IhdWKIBo0aBDx8fG0bNnSqD0lJQV/f3+jX1p+fn4kJydXeU1Wds54tnsJgLnvvMPx48erfJvmbNGiReTm5lLfpy11PAKqdFsqlQWeQYMNf0hIf4JCKKtWBJGHh8d92/Pz87G1Ne5/zNbWloKCguooCxfvILyaP4W2sJBRo0Zx8+bNatmuublw4QJr165FZWGBW4u+1bJNKztnGgQNBOD999/nt99+q5btCiFKqxVB9CD29vZoNBqjNo1Gg729fbXVUL91P+o1bM7NmzcZPXo0hYWF1bZtczFv3jxKSkrwbBaOnXODattu3UbtaBDQGa1Wy//7f/9POrcVQiG1Ooj8/f25ePGiUVtqaipNmzatthos1JY0eGIkdeq6c+zYMUaNGiV/XZfD999/z/fff4+1jS31mv+j2rfv0bofLh4+XLp0ifHjx1NSUlLtNQjxqKvVQRQaGoparSYmJgatVktCQgKJiYn07t27Wuuwsq1Dw9DRWNvYsXv3bnr27MngwYPZuHEjx48fly5lHuDSpUtMmjQJAM+WfbCyc672GiwsrfEKGWW4XzRv3rxqr0GIR12tDiJra2tWrlzJnj17CAsLY8GCBURHR+Pt7V3ttdi7ePNYz5l4tYjAytqWPXv2MGXKFPr06UNAQADR0dHVXlNNptFoePnll7l79y4evu1xa9ZVsVpsHNxoFBKFhYWa//3f/5VhxYWoZiq9vARj5OrVq0RERODb7Q2s7F0rtI7iwlwyfk9El5uG5u51sm5dQafTMXfuXEaNGmXiimu2gwcPkpycTHBwMIGBgahUKpKSkli+fDnx8fE4uXjg03UaltbVd1/vQW4nJ3D5l024urqye/du6tWrp3RJQjwSLJUuwBxZ2jhSP+C/PQNkXPyJ3xNXM3v2bNzd3enbt3qeDFPStWvXmD17ttHZhbOzMxYWFoaHAiwtrfB6YmSNCCEAt6ZPkX/zDLevnmHGjBnExMRU6ftMQoh7avWludrC1TcM73bPodfr+ec//8mPP/6odElVasuWLTz11FPs2LEDaxtbPJuG4ORSn+zsbDIzM6njXA+vxzrg2/nVSo81ZEoqlQUeQf+DtY0dO3bsYOvWrUqXJMQjQS7N/YUpLs3dj16v5+apzVw7m4BKpWL69OlMmDDB7P7i/umnn3jppZcoLi7GwzcY91bPY+1w7/tYmHcH9HqsHdxq9H7fSTnApSPrcXZ2Zs+ePdSvX1/pkoQwa3JGVE1UKhUebfrh3e5ZAD744APGjBlDXl6ewpWZTlpaGmPGjKG4uBivlj1oGDrKEEJw76EAG8d6NTqEAFz9OuLeuDXZ2dksX75c6XKEMHsVDiKdTseFCxdITEzk6NGj/P7779L550OoVBa4N++FX6ex2No5sGPHDvr27cvvv/+udGmVlp+fz6hRo7hz5w7u3q3waFl774OpVCrcWvQBYM2aNdy5c0fhioQwb+V+WOHo0aOsXr2agwcPGv01r1KpcHJyolOnTgwaNIj27dubtFBz4uzVGpuuU7n2cyznzp3jH08/zcIFC/Dx8aGwsBAXFxd8fX2VLvOhLl26xMqVKzl+/DhnzpxBq9Xi5OJBg8eHo7JQK11epdi7eOPh046bv58gJiaGGTNmKF2SEGarzPeILl26xKxZs7h27Rrdu3enc+fO+Pv7U7duXXQ6HZmZmZw7d44jR46wY8cOvL29mTt3bq34hfpnVXWP6H5Kigq4cWwdNy8eKzWtW7duTJgwodT4SzXF9u3bmTJlCjk5OYY2Vw8fPNoPwc7ZS8HKTCfvzkWSflyAg4MDiYmJuLi4KF2SEGapzEE0YMAAxo8fT5cuXR46r16vZ9euXXz66ads2rSp0kVWp+oMIgC9Xsed8/8hN+0XVCoLVGpLsm+nUaS914deREQEsbGxWFtbV3ktZVFQUMD777/P559/DoCHb3ucfDph79IYtbWdwtWZ3tVDy0m//CuTJ09mypQpSpcjhFkqcxDp9fpy32SuyDJKq+4gup/iwlzuXNjLnQt70BTkMWPGDF599VXDdJ1Oh0ajMQwS6OxcNV3jFBYWUlJSgr29PTqdjvj4eObPn8/169dRqy1pGNQPN/8ute4Yl0du+gXOJyzC0dGRfv360bZtW8LDw2nUqJHSpQlhNsr8sMKff9n861//IikpqdQ8GRkZRERE3HcZUXaWNo54tOxNo9Ao4N5YPZcuXQLg9OnThIWF0axZMwIDAwkMDGT69Okmf1Dkxx9/pG3btjRr1oxWrVoRFhbGP//5T65fv46rhw9Nu02mXtOnzP4YO9ZviqdvELm5uaxevZopU6bQqVMnNm/erHRpQpiNCj01Fx8fz8CBA/n222+N2nU6HdeuXTNJYQKcPFvg2TSMwsJCZsyYwS+//MKAAQNIS0vD0soaWzsHLCzUrFu3zmiU2G+++YaePXuydOnSUsNklMXGjRsZNWoUOTk5WFioyczMJC0tDQcnV3xCh9Hkqak4uPmYbkdruAZPRNGs62s0bPc89X3aotVqmThxIvPnzy81uqter+fy5cvS0a0Q5VChF1qbN2/OjBkzWLx4MQMGDOD111/HwsKC27dv06lTJ86ePVsVtVaLmnBp7s+KNDkk73wHTX4uVlZWFBUV4eH3OA2Ch2GhtiTryjFSD8ZiZWXFli1b+OWXX5g9e7ZheW9vb95880169+790LMXnU7Hxx9/zMKFCwFo1PofuLfoTXFhDkWau9jW8cDCsmbcq1LS7Qt7uXpsEzpdCcHBwXTv3p2QkBB+/fVXNm7cyNmzZ3F2dmb48OFERUXh7u6udMlC1GgVDqKDBw+Snp7O+PHjadSoEYsXL0av1xMeHi5BZGJ3Ug9y6ed1AHg2DaNB+8FGj0ffOLmJa2cTcHR0JDc39958zbtRmJFMZvoV4N6QGW+//TatWrXi4MGDxMbGkpmZSVRUFM888wwZGRm89tpr7N69G5VKRePHX8LNv3P172wtcff6Ga78/AWFBaVfSLaytjU8bGJjY0NERAQ9e/YkIiICV9ea8X9KiJqkQkHUokULDhw4gJubm+EX2OXLl5k3bx5RUVESRCam1+u5c247WFjh1qwHKgvjK6q6kiIu7/+YjBupWFioaRwyGFefMPS6EjIuHuLW2e3k595FpVLh4+NTajDBFi1acOfOHdLT07Gzr4PX40Nw9mpdnbtYKxUX5nL35jm0GakUZv2OpZ0z9l5P4OzVhvzMS2Sn7DZ6NN/S0pK3336bESNGKFe0EDVQpc6I3NzcACgpKWHevHl89dVXFBUVSRApQJufSeb5ndh5tsXJs7nRtGJtPneSdnLjXAIlJcXYOzrj6tcZC2sHMs7/SO7dez0H1GvYHM/goVjby/sypqLNyyD72mm0t89w89Jp4N7DJy+99JLClQlRc1RoGIhXX30Ve/v/dt2vVqt58803admyJfHx8SYrTpSdtb0LHu3u/8vN0toej9bPUdc3HM3dm9Sp38xwr8fVtwMZFw+j0+lwb9q51NmWqBxrB1fcm3WBZl2wqbeby7/EMXXqVNRqNW5ubhw4cID8/HzDJW4hHkVlPiO6cOFCmVfatGnTChektNp6RiRqh1vnvufKiW9KtdetW5elS5fSteu9kWpr4zt4QlRUmc+I+vTpg0qlMvyA/PkH5c9ZplKpavWlOSGqUr2AXuhLirj22w+4efph5dKUotxr3Lp0iqFDh/LMM89w69Ytzpw5g5ubG+vWraNx45ozZpMQVaHMZ0RpaWlGn/V6Pc888wwxMTF4eRn3LdawYUPTVVjN5IxIVAfjP+R03D63k6untpV6MdnLy4u4uDh8fHwUqFKI6lHmM6IHhYunp2etDh4hlPDny24qlQXuLSKxdw8g93Yqtk6e2NRx5+bxdVy7lswLL7xAXFwcfn5+ClYsRNUxyzvTxcXFTJkyhUGDBvHPf/6TwsJCpUsS4qEc6vni0TwCZ6+W2NapT6MOr+Dm9Rg3btygb9++fP/990qXKESVMMsg2rlzJx4eHqxfv542bdqwZcsWpUsSotzUVrY06jAW98atyczMZNSoUUybNs2sRvUVAsw0iE6cOMGTTz4JQMeOHfn5558VrkiIirkXRi/TOLg/llbWrF+/npCQEObOnUtqaqrS5QlhEmW+R7Ru3bpSbSUlJWzdurXUgGGDBw+ufGWVkJubi4ODAwAODg7yF6So1VQqC+o164qDewDpp77kzrVkVqxYwYoVK3jsscdo3bo1QUFB9O/fH0dHR6XLFaLcyhxEn332Wam2evXqlXqBVaVSKR5Ejo6Oht6P8/Ly5IdTmAW7ul406TyJend+5+6lg9z+/Sjnz5/n/PnzbN68mW+++YavvvoKKysrpUsVolzKHEQJCQlVWYdJtWnThkOHDtGxY0cOHjxIUFCQ0iUJYTIObj44uPng0XYABVlp5GdcJiN5Jz///DMffvghM2fOVLpEIcqlzEHUo0cPwsPD6dixI2FhYTX6LCMyMpI9e/YwcOBAXFxciI6OVrokIUzOQm1lCCU7Zy8u7PmY5cuX88QTT9CzZ0+jeTMzM/n555+xsbHB0dERX19fQ1+RQiitzC+0Xrx4kUOHDnHgwAF++eUX/Pz8CA8P58knn6Rdu3ZYmEkfZfJCq6itbiX9yJXjW3ByciYu7itatWoF3Ht4Z+TIKNLTbxrmtba2YfHiaJ599lmlyn2gQ4cO8d5779G1a1defvll6tSpo3RJoopVqPftkpISjh8/zqFDhzh48CC///47jz/+OB07dqRjx440adLE5IWeOnWK0aNHGz0Bl5SUxOzZszl37hyenp7MmDGDLl26VGo7EkSittLrdaQlxpL++wksLCwYPHgwbdu2ZebMmRQWFuJSvwk29k4UaXK4c+N3AP71r38xbty4B/Zrd/LkSd555x2eeOIJJk6ciJ2dXZXuw/HjxxkwYIDhHq+LiwsTJ05k5MiRcu/LjFUoiP4qJyeHw4cPc/DgQQ4fPszOnTtNUZvBjh07mDVrliEAAbRaLZGRkQwdOpQhQ4awb98+pk6dyrffflupnh4kiERtVlKk4daZb7mRtBedrsTQ3iCgMx5tXsRCbYler+d28m6uHt+MXq/n+eefZ/To0bRt29YokL788kvemDED7f+9EN6kSRPmzZtX6o+9o0ePsmrVKqysrHB2dqZRo0b07duX+vXrA3Dw4EEWLFjA7du3CQ0NpUOHDnTu3Nkw/Q/nz5/n+eefJysrC3efIPRF+dxOSwKgXbt2LF26FF9f3yr5vgllmSSIqlJ0dDT79++nT58+fPLJJ4Yg2r9/P6+//joHDx40/PC8/PLLBAUFMXbs2ApvT4JImIOC7OvcOfM1d9LO0bDNs7g161rqrCfr6nGuJP6boiItAP7+/oSEhGBtbc2tW7fYsWMHAB5NO6C9m0Zm+mUAxo8fz4wZM1CpVCQnJ9OnTx/DyMB/sLS0pFevXhQVFd33D1MLCws6d+5Mv379sLCw4OTJk2zdupX09HTq+7SjYcgoUFlw99qvpJ/eRE7WLezt7Xn33XcZMGCA9ExuZsr0sEJKSgrp6em0adPG8H4OwO7duw3d1leVQYMGMWnSJBITE0vV5O/vb/Qf0s/Pj+Tk5CqtR4jawM65AY06vEJDXYnRsPJ/VrdREDZ1PMn+/RBZV34hJSWFlJQUw3RLSysaBr+Em++T6HUlOCYncO3UNyxbtoyioiImTJjA8OHDyc3Nxd0nCLv6LSnR5lFy9zLpv59g+/btAFjb2FK/eU8c6geQdyuF4uwUbl3+jT179rBnzx6jmuo1aoHX4yMMNTs3bI2Duz/pp77ixoWfmTx5MqdPn+btt99Grb7/fona56FBtHr1atavX4+fnx9nz55lxowZhidylixZUuVB5OHhcd/2/Px8bG1tjdpsbW0pKCio0nqEqE0eFEJ/sHNugF3bF/Bo/Rw5N5PQ5meg15Wg15VQx7M5ds5ehvW4B/TAysGdSz99TkxMDF9//TXp6em4ePjQ8PHhhsEWAeq1yiLr0k+gL8HFNxwrO2cAHNx8ge7Ub5dL5pVf0N76DZXaCguHRti5NsbJo3mpmi2t7WkQPBxbtwCu/rKRL774gqtXr7J8+XKjATpF7fXQIIqLiyM+Ph57e3uuXr3KxIkTSUtLY+TIkaW6rK9O9vb2aDQaozaNRiP/MYWoAJWFGqcGgQ+dr26jdtBhNJcOxZKeno69Y10ahow2CiEAa/u61G8R+cD1WNo44t60CzQt28NFKpUKV98OWDu4c+Wnlfz4448899xzzJ07l7CwsDKtQ9RcDw0ivV5v+OXeqFEj1qxZw8SJE7lx44aiQeTv709sbKxRW2pqqry8KkQVq9uwDRbhY8m9fADnpj2xdqi+e6mO9Zvi+9Rkrhz+lN9++40XXniBJ598kqeeeooTJ05w7NgxHBwc6N+/PwMGDHjgFRVRszz05R9XV1ejEVcdHBxYsWIFmZmZnD9/vkqL+zuhoaGo1WpiYmLQarUkJCSQmJhI7969FatJiEeFU4NAvELH4ODmU+3btnXywK/rdBq16YONnQOHDh3i/fffZ8eOHdy4cYOUlBTmz5/PE088wZQpU8rU1+SlS5eIjo7m008/ZdOmTZw+fboa9kT84aFPzd24cQO1Wo27u3upab/88gvBwcFVVtyfJSYmMnbsWMNTc3Dvcc85c+Zw7tw5PDw8mD59eqXvWclTc0LUHsXafDJTD6DTZGHp5I1DPV8Kc2+TdzWR9Esn0JWU4Ofnx4oVKwgMvP+lx8TERKKiosjKyjJqf++99xgxYkTV74Qo3+PbW7durZFvYpuSBJEQ5qEg+zrXj64i69YVrG1smD1rFsOHDzd60nbz5s1MnToVrVaLW6OW2Dh5oivM5kbKUdRqNRs3bjQMKXPhwgUKCgpo3bp1herRaDRkZmYC0KBBg8rvoBkpVxC1bNmS4OBgZs+ejb+/f1XWpRgJIiHMh65Yy83T8VxP2gfAU089xcKFC0lLS2P58uX88MMPAHi16IZH636o/q+rspu/biXt1x9wdXVl1apVbNiwgY0bN6LX63nllVeYPn061tbWD9yuYfs6HStXrmTJkiVGZ1y9e/dmzpw5eHl5mX6na6FyBdG5c+d45513OHnyJMOHD+fVV1+t8i4/qpsEkRDmJ/PyMW6c+JKC/BxsbGwo/L/eIiwtrfBq+xz1mhlf0tfrdFz9aQW3Lv/3XpGFWg160OlKCAoK4qOPPuKxxx4znGFpNBrOnTuHo6MjTZo0IScnh9dee41du3YBoFZbYutQB60mnyJtIfb29kyZMoXRo0djaVnm/qfNUoV6Vti6dSsLFy5ErVYzY8YMevXqVRW1KUKCSAjzVFSQzY0TG7l16SS2dg64+XfC1b+L4R2nvyrW5nNp7yKy71zDw7c9bi2eoViby/Wjq8nJugXce5irbdu2ZGVl8euvv1JUVATc61nCxsaGvLw8bO0d8Wo/GOeGbVCpVGjzMrj12xZupv4CQFBQEIsXL6Zp06bV842ogSrcxU9eXh5Lly5lzZo1hIWFMWvWLBo3bmzq+qqdBJEQ5kuv15OfeRnbOvVRWz38ak5JUQFFBdnYOnka2ooL87j929dkpp2mIO+uoV2lUlG3XkNKirXkZN1Cr9fj5tWMBsHD7/uIe/a1X7l58ktys+9gbWPD/xs9Gi8vL+zs7AgMDKzwvajaqNJ9zaWkpPDee+9x9OhRRo8ezdixY8t07bSmkiASQpSFXq9Hm5dBfuZl1Fa2OLj6oLa+F266Yi1FmhysHVxQqR78lkyxNp/0X7dw4/zBUtPGjx/PtGnTTNbruE6nY/PmzRw7dgxXV1fc3d3RaDSkpKSQmpqKs7OzYQSFgICAh/bnd+nSJZYsWcLp06fp27cvI0aMqPA4dSbr9PTHH39k3rx5qNVqfvzxR1OsUhESREKI6nb3xlny08+BTou+uIAbqUfR63QEBwezfPlyGjVq9LfL6/V6Lly4QFJSEsHBwaWeyrt48SLTpk3j8OHDZarH29ub559/nn79+tG0aVNDKOXm5nLs2DG++eYb4uLiKC4uNixTt25dfvvtt3Lu+T0VDqK0tDTOnz/P+fPnSUpK4vz581y8eBGdTmf0AmxtI0EkhFBa7q0LXDv6b3Kz7+Dk7Ez0okVERhp3mVRQUMDu3bv57rvvOHDgAOnp6cC9S4RhYWF06dKFrKwsLl++zK5duygsLMTe0Zm6fl3Q64pRFeeBygILO3ds69SnqCAbbdYF7l4/R15OpmE7tra2uLq6Ym9vT2pqKjqd7t52LCzwbNoB23otyLm0j9tp50lLS6vQ/pYriGbNmsX58+dJTk42DFzl7e1N8+bNCQgIoHnz5jRv3rxS4wEpTYJICFETFBfmcuP4etJ/PwFAVFQUffr04fTp0xw7doxdu3YZDb/hUKcuDi5eZFxLpri4qNT6PJt2oH7r57G0+fvLZ3qdjtxbyeSm/ULmlRNoCv67DQu1GjcPX6zq+lLXtyO2de6NKaXX68m5eY69696o0L6WK4heeuklo8AJCAgwGhbCHEgQCSFqCr1ez50Le0g7sYWSkuJS090a+GHn0ZY6DVph6+SJSqWiRFtAVtpJiu5eAas6WNu7YVu3AfZ1//7y3oO2ryvRUqzJpaQoH5s69VFb2jxw/m0fVazDg3I9vP7ll19WaCNCCCHKT6VSUa9ZV+zc/Mg4swVdSRE2To1QO3pRx7OF4Yzkz9TWdrj5hgGV75VcpVKhtrRB7WgDuFV6fQ9Sobeovv76a1q2bIm/vz8WFg/tN1UIIUQlOLg2wSH8NaXLqDIVCqL58+eTnZ2NtbU1zZo1IzAwkBYtWhAYGEhAQECpAeuEEEKIB6lQEP30009cv36d3377jbNnz3Ly5Eni4+MpKSnBwsKiwo/wCSGEePRUuIOjBg0a0KBBA7p37w7cGy5i0qRJZt87txBCCNMy2Q0eT09P5s2bx969e021SiGEEI+ACgVRRkbGfdt9fHw4c+ZMpQoSQgjxaKnQpbknn3wSDw8PwwMKgYGB+Pn5cerUKUPvs0IIIURZVCiINmzYwNmzZzl79ix79uwhNjYWrVaLWq1m+vTppq5RCCGEGatQEAUFBREUFGT4XFJSQlpaGs7Ozjg7339sDyGEEOJ+ynyP6O7duw+cplarady4cakQ+rtlqtOaNWuIiYlRugwhhBD3UeYgGjJkCDExMeTk5Dx03oyMDJYvX87gwYMrVZwpzJkzh9WrVytdhhBCiAco86W5DRs2sHjxYjp37kz79u3p1KkTzZo1w8XFBZ1OR2ZmJufOnePIkSMcOXKEfv36sWHDhqqsvUzCw8Np06YNt2/fVroUIYQQ91HmIHJwcGDmzJmMGTOGr776iu3bt3PmzBlKSkrurcjSksDAQLp06cI777yDh4dHlRX9V3FxccTFxRk+u7u7s2zZMgC6d+9OfHx8tdUihBCifMr9sIK7uzvjx49n/Pjx6HQ6srKyUKlUuLi4VEV9ZdK/f3/69++v2PaFEEJUXIW7+AGwsLDA1VXG7BFCCFFxFQqi5s2bG8Yw/ysrKys8PDx4+umnmThxImq1ulIFCiGEMG8VCqK5c+fy8ccfM2HCBNq1awfA6dOnWbJkCQMGDKBZs2YsXboUgEmTJpms2Irq16+f0iUIIYR4gAr1NffZZ5/x3nvvMXDgQMOw4f379+f9999n+/btREZG8vbbb/P1119XuLBTp04REhJi1JaUlMTAgQNp164dkZGR0sGqEEKYgQoFUXp6Ot7e3qXaPT09uX79uuHf2dnZFSpqx44dREVFGfVbp9VqeeWVV+jVqxdHjhxh2rRpvPbaa6SlpVVoG0IIIWqGCgVR+/btWbBgAVlZWYa2rKwsPvroI9q3bw/Ajz/+iI+PT7nXHR0dTWxsLOPGjTNqT0xMRKPRMGLECKysrIiIiCAkJIRt27ZVZBeEEELUEBW6R/TOO+8wZswYOnfuTMOGDdHpdFy/fh0/Pz+WLl3Kvn37WLRoEUuWLCn3ugcNGsSkSZNITEw0ak9JScHf39/oIQk/Pz+Sk5MrsgtCCCFqiAoFkZeXF1u3buWnn34iKSkJS0tLmjVrRocOHQCwt7dnz549FXq0+0Evwubn52Nra2vUZmtrS0FBQfl3QAghRI1R4feI1Go1FhYWWFpaotPpKCoqori4GEtLyyp5t8je3h6NRmPUptFosLe3N/m2hBBCVJ8KBdGtW7cYP348Z86cMbo05+vry6pVq3BzczN1nfj7+xMbG2vUlpqaajQchRBCiNqnQg8rvPfee6jVanbt2sUPP/zAjz/+yH/+8x+cnJyYN2+eqWsEIDQ0FLVaTUxMDFqtloSEBBITE+ndu3eVbE8IIUT1qFAQ7d+/nzfffNPofo6npydvvPEG+/btM1lxf2Ztbc3KlSvZs2cPYWFhLFiwgOjo6Ps+Ri6EEKL2qNClORsbGywsSmeYhYWFoTfuygoNDeX48eNGbY899hjr1683yfqFEELUDBU6I+rUqRPvv/++0Rg/t27dYv78+YSHh5usOCGEEOavQmdE06dPZ/jw4XTt2hUvLy8Arl27RkBAADNnzjRpgUIIIcxbhYLIzc2NrVu3sm/fPlJSUrCxscHf358nn3zS1PUJIYQwc2UOogsXLpRq8/b2NnpY4I95mjZtaoLShBBCPArKHER9+vRBpVKh1+sNbX8dk0iv16NSqTh79qzpKhRCCGHWyhxEu3btqso6hBBCPKLKHEQNGzasyjqEEEI8oir0+LYQQghhKhJEQgghFCVBJIQQQlHlfo8oIyODPXv2kJycTF5eHo6OjrRo0YIuXbrg5ORUFTUKIYQwY+UKotjYWJYuXYpKpaJhw4Y4OTmRm5vLmjVrsLS0ZNKkSQwbNqyqahVCCGGGyhxEcXFxLF26lGnTptGvXz/s7OwM0zQaDfHx8SxcuBAPDw969epVJcUKIYQwP2UOorVr1zJ16lQGDx5capqtrS2DBg2ipKSE1atXSxAJIYQoszI/rHDp0qWH9qwdHh5OSkpKpYsSQgjx6ChzEGk0GhwdHf92njp16pCdnV3pooQQQjw6yvX49l/7lhNCCCEqq1xPzcXHx2Nvb//A6Xl5eZUuSAghxKOlzEHk5eXFhg0bHjpfgwYNKlWQEEKIR0uZgyghIaEq6xBCCPGIqtAIrbVBdnY2U6ZMobCwEGtraxYuXIiLi4vSZQkhhPiLcj2skJuby9dff01ubi5wbyC8lStXMmzYMF577TVOnTpVJUVWxNq1a+nVqxdr1qwhMjKStWvXKl2SEEKI+yhzEF2+fJl//OMfzJ07l8zMTADee+89Fi1aRL169XBycmL48OEcO3asyootj2HDhtG3b18AioqKsLGxUbgiIYQQ91PmS3Mff/wxLVu2ZPHixdja2nLnzh02bNhAZGQkixYtAqBJkyYsXbqUzz//vMoK/qu4uDji4uIMn93d3Vm2bBl16tQB4Ny5c2zYsEHOiIQQooYqcxAdPnyYFStWYGtrC8D+/fvR6XQ899xzhnk6duzI8uXLTV7k3+nfvz/9+/e/77Rjx44xe/ZslixZgrOzc7XWJYQQomzKHEQ5OTm4ubkZPicmJmJpaUlISIihzcHBAZ1OZ9oKKyg5OZk5c+awYsUKvLy8lC5HCCHEA5Q5iBo2bEhqaipeXl6UlJSwf/9+nnjiCaNeuH/66ScaNWpUJYWW1yeffIJGo+H1118HoFu3bowcOVLhqoQQQvxVmYOoX79+vPPOO7z66qscPnyY27dvM3v2bMP0o0ePsnjxYoYOHVolhZbXkiVLlC5BCCFEGZT5qblRo0YRERHBu+++y+7du5k8eTI9evQA4N1332XIkCEEBwczatSoShd16tQpo0t+AElJSQwcOJB27doRGRnJ3r17K70dIYQQylPp9Xp9ZVeSlJSETqejRYsWlS5ox44dzJo1i5KSEo4fPw6AVqslMjKSoUOHMmTIEPbt28fUqVP59ttvadiwYaW3+WdXr14lIiIC325vYGXvatJ1CyGEOdv20bMVWq7MZ0T5+fm89dZbhISEEB4ezpw5cwwvtgYEBJgkhKKjo4mNjWXcuHFG7YmJiWg0GkaMGIGVlRURERGEhISwbdu2Sm9TCCGEssp8j+iTTz5h9+7djB49GrVazbp168jKymLx4sUmK2bQoEFMmjSJxMREo/aUlBT8/f2NhqHw8/MjOTnZZNsWQgihjDIH0XfffceCBQvo0KEDAE888QSDBw+mqKgIKysrkxTj4eFx3/b8/HzD+0t/sLW1paCgwCTbFUIIoZwyX5q7desW/v7+hs+tW7dGr9dz586dKinsz+zt7dFoNEZtGo3mb8dGEkIIUTuUOYhKSkpQq9WGzyqVCisrK4qKiqqksD/z9/fn4sWLRm2pqak0bdq0yrcthBCiapWr922lhIaGolariYmJQavVkpCQQGJiIr1791a6NCGEEJVUqaHCS0pK2Lp1a6lxfgYPHmya6v6PtbU1K1euZM6cOXz66ad4eHgQHR2Nt7e3SbcjhBCi+pX5PaJu3bqVbYUqFbt27apUUUqS94iEEKJiKvoekQwVLoQQQlG14h6REEII8yVBJIQQQlESREIIIRQlQSSEEEJREkRCCCEUJUEkhBBCURJEQgghFCVBJIQQQlESREIIIRQlQSSEEEJREkRCCCEUJUEkhBBCURJEQgghFCVBJIQQQlESREIIIRQlQSSEEEJREkRCCCEUJUEkhBBCURJEQgghFGWpdAFVpbi4mOnTp3Pjxg08PDz44IMPsLa2VrosIYQQf2G2Z0QJCQl4eHiwfv16/P392bp1q9IlCSGEuA+zPSPq2bMnERER6PV60tPTadmypdIlCSGEuI9aH0RxcXHExcUZPru7u7Ns2TIA1Go1Q4YMIT09nVdffVWpEoUQQvwNlV6v1ytdRFXbs2cP27Zt46OPPnrovFevXiUiIgLfbm9gZe9aDdUJIYR52PbRsxVazmzvEX399desX78eADs7O9RqtcIVCSGEuJ9af2nuQbp3787UqVP57rvvsLKy4p133lG6JCGEEPdhtkHk6OjIp59+qnQZQgghHqJGXpo7deoUISEhRm1JSUkMHDiQdu3aERkZyd69exWqTgghhCnVuCDasWMHUVFRFBUVGdq0Wi2vvPIKvXr14siRI0ybNo3XXnuNtLQ0BSsVQghhCjUqiKKjo4mNjWXcuHFG7YmJiWg0GkaMGIGVlRURERGEhISwbds2hSoVQghhKjUqiAYNGkR8fHypl09TUlLw9/dHpVIZ2vz8/EhOTq7uEoUQQphYjQoiDw+P+7bn5+dja2tr1GZra0tBQUF1lCWEEKIK1aggehB7e3s0Go1Rm0ajwd7eXqGKhBBCmEqtCCJ/f38uXrxo1JaamkrTpk0VqkgIIYSp1IogCg0NRa1WExMTg1arJSEhgcTERHr37q10aUIIISqpVgSRtbU1K1euZM+ePYSFhbFgwQKio6Px9vZWujQhhBCVVCN7VggNDeX48eNGbY899pih7zghhBDmo1acEQkhhDBfEkRCCCEUJUEkhBBCURJEQgghFCVBJIQQQlESREIIIRQlQSSEEEJREkRCCCEUJUEkhBBCURJEQgghFCVBJIQQQlESREIIIRQlQSSEEEJREkRCCCEUVSOHgVBSSUkJAEUFWcoWIoQQtczVq1fx9PTE0rJ80aLS6/X6KqqpVjp69CiDBw9WugwhhKiVdu3aRaNGjcq1jATRX2g0Gn799Vfc3d1Rq9VKlyOEELWKnBEJIYSodeRhBSGEEIqSIBJCCKEoCSIhhBCKkiASQgihKAkiIYQQipIgEkIIoSgJIiGEEIqSIDITX331FS1btiQoKMjwtWXLFqXLMolTp04REhJi+KzVapk9ezahoaGEhoYyb948Q9dMtdlf91Oj0ZQ6plFRUQpWWDmnTp1i0KBBPP744zz11FN88skn6PV6szueD9pPczqeu3fv5plnniEoKIgePXqwceNGoBI/m3phFt566y39okWLlC7D5LZv364PDg7Wt2vXztC2cOFC/aBBg/SZmZn669ev65999ln9ihUrFKyy8u63n8ePH9d36tRJwapMJzc3Vx8WFqZfs2aNvri4WH/x4kV9t27d9GvXrjWr4/l3+2kuxzMtLU3fqlUr/f79+/V6vV5/9uxZfevWrfUnT56s8LGUMyIz8dtvv9GiRQulyzCp6OhoYmNjGTdunFH7li1bGDNmDHXr1sXT05OxY8eyefNmhaqsvAftpzkd0+vXr9O+fXuGDBmCWq3Gx8eHHj16cOzYMbM6nn+3n+ZyPL28vDh8+DDh4eHodDoyMzNRq9U4OjpW+FhKEJmBoqIizp8/T3x8POHh4fTo0YOYmBj0tbz3pkGDBhEfH0/Lli0NbXfv3uXWrVv4+/sb2vz8/Lh06RJarVaJMivtfvsJ94Lo1q1bPPPMMzz55JNMnDiRmzdvKlRl5TRt2pRly5YZPmu1Wvbt20fz5s3N6nj+3X6a0/F0dHQkNzeXVq1aMWLECIYOHUq9evUqfCwliMxARkYGbdq0oV+/fiQkJPDxxx+zfv161q5dq3RpleLh4VGqLT8/HwA7OztDm52dneEafG10v/0EsLe3Jzg4mH//+99899132NjY8Morr1Rzdaan1WqZPHky1tbWPP3004B5Hc8//Hk/Bw8ebHbH087OjhMnThAXF8eXX35p+H1TkWMp4xGZAQ8PD9atW2f4HBgYyNChQ9m5cydDhw5VsDLT++M/+Z//YxcUFAD3fnGbkzfffNPo84wZM+jQoQNXr14tdzf7NcWtW7eYMGECAF988QUWFvf+Fja34/nX/bS3tze746lWq1Gr1bRp04YXX3yRU6dOARU7lnJGZAaSkpJYunSpUVthYSHW1tYKVVR1nJ2dcXd3JzU11dCWmpqKj49Pubuer8n0ej3R0dFG+/nH5Q0bGxulyqqU8+fP88ILL+Dj48Pq1atxcXExy+N5v/00p+N5+PBhBgwYYNRWVFSEk5NThY+lBJEZcHBwICYmhs2bN6PT6Th16hRr167lxRdfVLq0KtG3b1+WLVvGnTt3uHnzJp9++inPPfec0mWZlEql4rfffuPDDz8kJyeH7Oxs3n33Xbp06YK7u7vS5ZVbZmYmUVFR9OnTh/nz5xv9kWROx/NB+2lOx7NFixZcvnyZ1atXU1JSwtGjR9myZQsvvvhihY+ljEdkJg4cOMCiRYu4ePEiLi4ujBo1ymxGmk1MTGTs2LEcP34cuHe298EHH/DDDz9QXFzMs88+y+uvv17rBzL8637euXOHd999l0OHDlFSUkKXLl2YNWsWzs7OCldafl988QXz58/Hzs4OlUplaO/UqRMLFiwwm+P5d/s5e/Zsszmep0+f5r333iM5ORkvLy8mTJhAz549K/yzKUEkhBBCUXJpTgghhKIkiIQQQihKgkgIIYSiJIiEEEIoSoJICCGEoiSIhBBCKEqCSDxyunXrRkBAwH2/HnWJiYmG78XcuXO5evUqAQEBnD9/vtS88fHxhIaGPnSd0dHRhnXu3r27KsoWtVzt7ENDiEqaPHky/fr1U7qMGus///kPLi4uZGVlVXpdY8aMYciQIYSHh1e+MGGWJIjEI8nBwaHWda1SnVxdXXFwcDBJEDk4OODg4FD5ooTZkktzQvzFJ598wujRoxk1ahTBwcGGIdc/++wzunbtSlBQEP/zP//DiRMnDMsUFxezcOFCwsLCCA0NJSYmhn79+hEfHw/A0KFD+eCDD4y289dLVX+3/jfeeIM5c+bwxhtvEBQURIcOHYzGvdHpdCxbtoynnnqKoKAghg0bRkpKCunp6QQGBnL48GHDvHq9nm7durFp0yaTfc/i4+Pve6nT3Hp/F1VDgkiI+9i/fz8hISF89dVXdOnShY0bN7J27VrmzJnDli1b6NKlC8OHD+fq1avAvfDasmULH374IatWrSIhIYHk5OQyb+9h6wfYtGkTXl5ebN26lZEjR7JkyRLOnDkDwNKlS1mzZg3/+te/2LJlC+7u7owZMwY3NzfCwsLYsWOHYT3Hjx/n9u3bREZGmui7BU8//TQHDhwwfC1evBhLS0tGjx5tsm0I8yVBJB5J8+fPJygoyOjr6NGjhul2dnaMGTMGf39/XF1dWbFiBZMnT6ZLly74+PgwduxYgoODWbduHXq9nq+++opx48bRuXNnWrRowUcffURJSUmZ6/m79f+hcePGTJw4kcaNGxuGYz59+jR6vZ4NGzYwbtw4evbsiY+PD7Nnz6ZHjx5kZ2fTt29fdu7cSVFREQDffvstXbt2xdHRscz19e/fv9T3a/bs2Ybptra2uLu74+7uTnFxMe+88w7jx4+nS5cuZd6GeHTJPSLxSHr55Zfp27evUdufR0pt2LChoffkvLw8rl27xptvvsmsWbMM82i1WqytrcnIyCAjI4PAwECj5T09PctUy8PW/4fGjRsbLefg4EBxcTGZmZlkZGTQunVrwzQnJyfeeOMNAHr06MGcOXM4dOgQ4eHhfP/998ydO7dMtf1hyZIl+Pj4GLXt3LmT2NhYo7bCwkLGjx9Pu3btavXoo6J6SRCJR5KLiwtNmjR54PQ/D1b2x5nN/PnzjcIG7p0J2NraAvfu0/zZ3w1MWFxcXOb1/9369Ho9VlZWhn/fj4ODAxEREXz//fdYWloahiAojwYNGpT6frm5uZWa76233iI/P58PP/zQaBgEIf6OXJoT4iH+GHny5s2bNGnSxPD173//m/379+Pg4EDDhg05efKkYZnMzEzS0tIMn62trcnLyzN8vnLlSpnX/zB16tTBzc3NcL8I7g3X3LFjR8PwzX379mX37t3s2rWLXr16GcLLlFatWsV//vMfli1bVq7LfkLIGZEQZTB69GiWL19O/fr1adWqFd988w1ffvkla9asAe69K7No0SK8vb3x9fVl4cKFhqGgAVq1akVcXBzPPvssDg4OpUYpfdj6H2bEiBH87//+L40aNaJJkyYsX74cR0dHmjdvDkB4eDhqtZpNmzbx+eefm/A7c8/hw4dZuHAh77//PnXr1uXWrVuGafKYvHgYCSIhymDYsGFoNBoWLFjA7du38fX1ZcmSJbRv3x6AgQMHkpWVxZw5cygsLGTUqFGGkVYBoqKiSE1NJSoqChcXFyZOnGh0xvSw9T9MVFQUeXl5vPnmm+Tl5dG+fXtWrFhhCDu1Ws0//vEPEhISCA4ONuF35p6tW7dSVFTEtGnTSk1LSkoy+faEeZERWoWoIqGhobz++us1pgeHiRMn4uvry6RJkx44T2JiIsOGDePYsWMmfwk1ICCATz/9lK5du5p0vaL2kzMiIczczz//zJkzZ9i7d+99z1juJyMjA71eb5J7PXl5eeTn51d6PcJ8ycMKQpi57du3s2TJEqZNm4a3t3eZlunevTuLFi0yyfZjYmKknznxt+TSnBBCCEXJGZEQQghFSRAJIYRQlASREEIIRUkQCSGEUJQEkRBCCEVJEAkhhFDU/wcnLAIoo51XVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.signal import welch\n",
    "\n",
    "win = int(4 * sf)  # Window size is set to 4 seconds\n",
    "freqs, psd = welch(data, sf, nperseg=win)  # Works with single or multi-channel data\n",
    "\n",
    "print(freqs.shape, psd.shape)  # psd has shape (n_channels, n_frequencies)\n",
    "\n",
    "# Plot\n",
    "plt.plot(freqs, psd[1, :], 'k', lw=2)\n",
    "plt.fill_between(freqs, psd[1, :], cmap='Spectral')\n",
    "plt.xlim(1, 30)\n",
    "plt.yscale('log')\n",
    "sns.despine()\n",
    "plt.title(chan[1])\n",
    "plt.xlabel('Frequency [Hz]')\n",
    "plt.ylabel('PSD log($uV^2$/Hz)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen on the plot above, there is a strong alpha oscillatory peak around 10-11 Hz.\n",
    "\n",
    "## Extract the aperiodic and periodic component of the spectrum\n",
    "\n",
    "The Irregular-Resampling Auto-Spectral Analysis (IRASA) method ([Wen & Liu 2016](https://doi.org/10.1007/s10548-015-0448-0)) allows to separate the\n",
    "fractal and oscillatory components in the power spectrum of EEG signals.\n",
    "\n",
    "The steps are:\n",
    "\n",
    "1. Compute the original power spectral density (PSD).\n",
    "2. Resample the EEG data by multiple non-integer factors and their\n",
    "reciprocals ($h$ and $1/h$).\n",
    "3. For every pair of resampled signals, calculate the PSD and take the\n",
    "geometric mean of both. In the resulting PSD, the power associated with the\n",
    "oscillatory component is redistributed away from its original\n",
    "(fundamental and harmonic) frequencies by a frequency offset that varies\n",
    "with the resampling factor, whereas the power solely attributed to the\n",
    "fractal component remains the same power-law statistical distribution\n",
    "independent of the resampling factor.\n",
    "4. It follows that taking the median of the PSD of the variously\n",
    "resampled signals can extract the power spectrum of the fractal component,\n",
    "and the difference between the original power spectrum and the extracted\n",
    "fractal spectrum offers an approximate estimate of the power spectrum of\n",
    "the oscillatory component.\n",
    "\n",
    "**For more details on the IRASA method, please refer to**\n",
    "\n",
    "- Wen, H., & Liu, Z. (2016). Separating Fractal and Oscillatory Components in the Power Spectrum of Neurophysiological Signal. Brain Topography, 29(1), 13–26. https://doi.org/10.1007/s10548-015-0448-0\n",
    "- https://github.com/fieldtrip/fieldtrip/blob/master/specest/\n",
    "\n",
    "**For another method to extract fractal noise, please refer to**\n",
    "- https://github.com/fooof-tools/fooof\n",
    "- https://www.biorxiv.org/content/10.1101/299859v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02-Jun-22 14:35:06 | WARNING | The evaluated frequency range ends after the resampled Nyquist frequency (52.63Hz). Decrease the upper band (30.00Hz) or decrease the maximum value of the hset (1.90).\n"
     ]
    }
   ],
   "source": [
    "# Apply the IRASA technique\n",
    "freqs, psd_aperiodic, psd_osc = yasa.irasa(data, sf, ch_names=chan, band=(1, 30), win_sec=4, return_fit=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEgCAYAAABvi64JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJxUlEQVR4nO3deVxN2/sH8M9pkkqESBkiMkQaKENkSAMpMqXJkCEy5rpCpsslckUylHmWoTLE1W2QKZExdCmFMpYKTRrO+v3h6/ycW6icoU7P+/U6L87ae6/97LPrPO21116LwxhjIIQQQkRAStwBEEIIqT0o6RBCCBEZSjqEEEJEhpIOIYQQkaGkQwghRGQo6RBCCBEZSjoS7ObNm2jfvj2WLl0q9H1t3rwZdnZ2Vd5+7dq1cHZ2BgDExcWhffv2yMvLE1R4tVJeXh6OHz8ukLp+9fyKSmZmJlatWoUBAwZAV1cXFhYW2LZtG4qKinjrDBgwAO3bty/3tXnz5p/uIzQ0FO3bt0dgYGC5y0+cOAErKyvo6enB1tYWkZGRAjs+SUBJR4KdPn0ampqaOHfuHAoLC4W6r4kTJ2LXrl0CqUtfXx9XrlyBgoKCQOqrrfbs2YMjR46IOwyRSU9Px/Dhw/HixQusXbsWYWFh8PDwwNGjRzF37lzeeidOnMCVK1f4XpMnT4aMjAz69Onz0/18/b06efJkmWURERFYsWIFpkyZgtOnT2PYsGGYOXMm7t+/L9Bjrcko6UiooqIiXLhwAW5ubvj8+TMuXLgg1P0pKipCRUVFIHXJyclBVVUVHA5HIPXVVrXtue/ly5ejdevW2LZtG7p3744WLVrAwsICmzZtQkREBC5fvgwAaNiwIVRVVXmvZ8+eYffu3XB3d4eent4P9/Hu3Ttcv34ds2bNwrNnzxAfH8+3/MSJExg+fDiGDx+Oli1bYsKECTA2NsaZM2eEddg1DiUdCXXx4kV8+vQJpqam6NGjB99fZcHBwbC1tUVAQACMjIzQo0cPrF+/HqWlpbx17t27B3t7e3Tp0gXm5ubYsWMHuFwub/vhw4fjt99+g4GBAbZt21am+eXhw4cYN24cDAwMYGJiAh8fHxQXF/OWX7p0CUOHDoWuri7c3Nzw6dMn3rL/Nq+9evUK06dPh4GBAXr16oXVq1ejpKSk3OPOycmBp6cnjIyMYGRkhAULFiA3NxfAl0Ts7++PAQMGoEuXLhg7dizu3bvH29bZ2RmBgYFwc3ODrq4urKyskJCQgB07dqBnz54wNjaGn58f3/q+vr5wdXWFrq4uhg4dyvtiA7586R88eBAWFhbo0qULbG1tERMTw1vu6emJ5cuXw9PTE/r6+ujZsye2bNnCdzy7du1C//79oa+vj7Fjx+Lu3bsV2j44OBj+/v54+PAh2rdvj/T09HI/r29xuVxs2bIF/fr1g76+PlxcXPD06VO+5evWrYOxsTEMDQ2xdOlSvnO6d+9eWFhYoHPnzjAyMsLvv/+O/Px8Xjx2dnbYsWMHevfuDT09PcybN4+3/L88PT2/2wRW3rG8ffsWV65cgaurK6SlpfmW6enpYf/+/TAwMCiz3YcPHzB//nwYGhrCzc3tp59RWFgYFBUVYWlpCS0trTJXO9OnT8fEiRP5yjgcDt/Pd63HiERyd3dnDg4OjDHGjh07xtq3b89evHjBGGPs5MmTTEdHh40ZM4Y9evSIRUVFMSMjI+bn58cYYywzM5MZGhqybdu2sWfPnrGLFy+yfv36sYCAAN722trabMWKFezZs2fs5cuXzM/Pjw0fPpwxxlhqairT09NjK1asYMnJySwqKor17t2beXt7M8YYS0lJYTo6Omzz5s0sJSWFBQQEsPbt2zMnJyfGGGPXr19n2traLDc3l33+/JmZm5uziRMnskePHrGbN28yU1NTtnHjxnKP28nJiQ0fPpzdvn2bPXz4kNnY2LD58+czxhhbsmQJ69u3L7t48SJLTk5mS5YsYfr6+uzt27e8bbt27cpCQkJYamoqGzt2LOvevTubPXs2S05OZnv27GHa2tosMTGRt/7X40hOTmY+Pj5MR0eHpaSkMMYY2759O+vWrRs7e/YsS0lJYX5+fqxjx4687RcsWMB0dHTYpk2b2PPnz1lAQADT1tZmDx8+ZIwxduTIEdavXz928eJFlpqayrZt28Z0dXVZWlraT7cvKChg3t7ezMbGhr17946VlJT89Gdm06ZNzNjYmF24cIGlpqYyDw8PNmDAAFZSUsL8/PyYtrY2W7RoEUtNTWURERFMR0eHBQUFMcYYO336NDMwMGBRUVEsPT2dRUREMH19fbZnzx6+nzk3NzeWlJTELl68yHR1dXnL/+vjx4/s3bt35b7KO5aLFy8ybW1tlpWV9dPj/NbMmTOZkZERe/36dYXWHzZsGO/nacOGDUxPT4/l5uZ+d/3ExETWsWNHFhoaWqm4JBklHQn04cMH1rlzZ94v9Pv371mnTp3Ypk2bGGNfvgA6dOjA+/JijLE9e/awXr16MS6XyzZt2sQmTJjAV+fp06eZsbExb3ttbW2WkZHBW/5t0vH29mbW1taMy+XyloeFhTEdHR2Wl5fH1q1bx+zs7PjqnzhxYrlJJzo6muno6LD379/z1r148SI7fPhwmeN+8uQJX1JgjLF79+6xbdu2sQ8fPrCOHTuy8+fP85aVlpaywYMHsw0bNjDGviSRqVOn8pYfOnSIaWtrsw8fPvDK9PT02JkzZ3jru7i48MVgY2PD1q1bx7hcLjM2NmY7d+7kW+7q6so8PDwYY1+ShpWVFd9yIyMjdvToUcYYY/369WOnT5/mWz5hwgRe8v7Z9t+ek5/hcrmsR48ebN++fbyyDx8+sDVr1rD3798zPz8/1r17d1ZcXMxbPm7cOLZs2TLGGGOxsbEsPDycr86pU6eyhQsXMsbK/5lxd3fnfRa/6vTp00xbW5svvp85duwY09bWZhcuXKjQ+snJyUxbW5v9888/jDHGHj16xLS1tdmJEyfKXf/t27fMzMyMjR07tkJJv7aQEfeVFhG88+fPo6ioCObm5gC+tGF3794doaGhmDFjBgCgWbNmaN68OW8bXV1dZGZmIjs7G8nJyYiLi4O+vj5vOZfLRWFhIbKzswEACgoKaNy4cbn7T05ORteuXfnuyRgaGqK4uBjPnz9HUlISdHR0+LbR1dUt0z7+tS4NDQ00bNiQV2Zqavrd/crJyaF9+/Z89erq6uL+/fsoLS3lOyYpKSno6+sjKSmJV9aiRQve/+Xl5VG/fn0oKyvzyuTk5Ph6QnXr1q3McSQlJSErKwvZ2dll7hEYGhri77//5r1v2bIl33JFRUWUlJQgLy8Pr169gpeXF1/vw6KiIsjJyf10+8rKzs5GVlYWunTpwitTVlaGp6cn7726ujpkZGT4ln/+/BkA0KNHDzx8+BAbN25ESkoKkpOTkZKSgmHDhvHF9u3PjJKS0neb15YuXfrd+yBhYWFQV1fnK/t6P/Hjx498Pyvfk5qaitWrV2PMmDG835Ovtm/fjoCAAN77FStWwMbGBqdOnYKCggJMTEwAAB07dkSrVq1w8uRJjBgxgq+OtLQ0uLq6QlZWFv7+/mWa/GozSjoS6PTp0wAAMzMzXhmXywVjDLGxsQBQ5pfg6/0aKSkplJSUwNzcHHPmzClTd7169QCA74vvv+rUqVOm7Gv9paWl4HA4ZW5yy8rKllvX98oru+734uVyubzYAPB9qQL4aWeG/36OpaWlkJKSKvczKG9/5cXFGOPdX/P29kanTp34lsvLy/90+8r6+tn9aFspqbK3gL+uHxwcjOXLl8POzg59+vTBtGnTynQ/rsy5nD17NlxdXctd1qRJkzJlOjo6kJKSwv3799GvX78yy+fNmwdTU1PY2NiguLgY8+bNQ7NmzbBw4cIy69rb28PKyor3vlGjRmCM4ezZs8jPz+e7N8TlcvH8+XM8e/YMmpqaAL788TN+/Hg0atQIe/bsqVASrE2oI4GEefnyJW7duoWZM2ciNDSU9woODoaioiLvxufr16+RlZXF2+7evXto1qwZGjRoAC0tLaSmpqJVq1a8V1JSEjZv3lzuF89/aWlp4d69e3xfYHfu3IGsrCxatmwJbW3tMl1IHz16VG5dmpqaePnyJe8KCwBCQkIwatSoMuu2bt0aRUVFfFcu169fR//+/dG8eXPIysri9u3bvGWMMdy9exdt2rT56TF9z8OHD/nqS0hIQIcOHaCkpIQmTZrgzp07fOvfuXOnQvtTVlaGqqoq3r59y3ce9u3bx9dZ4Ucq0/uvXr16aNSoEd95KCwsRO/evSvU3XfPnj1wdXXF8uXLMWrUKHTo0AHPnz+vcg+6Ro0a8R33t6///mEAfLnSMTU1xa5du/iSOgDcunULZ8+eRf369QEAvr6+SEpKwoYNG1C3bt0ydTVo0IBvf0pKSoiPj8fLly+xcuVKvt+rvXv3QkpKivd79ebNG4wfPx5qamrYv38/JZxyUNKRMGfOnEGdOnXg4uICbW1t3qtTp04YPnw4IiIi8PHjRxQXF2PhwoVISkpCREQEAgIC4OLiAgBwdHTEs2fPsGrVKqSkpODKlStYtmwZ6tWrV6Gk4+DggPT0dKxatQpPnz5FTEwMvL29MWzYMCgrK2PMmDF4/vw51q1bh9TUVBw8eBAXL14sty4TExO0atUKCxcuxJMnT3Dz5k1s3rwZffv2LbOulpYWTExM4OXlhYSEBDx48ABr166FsbExlJSU4OTkBG9vb8TExODp06dYsWIFXr58idGjR1f5846KisLBgweRmpqKdevWIT09nVfflClTEBAQgHPnzuHZs2fw9/fH1atXeQ/B/sykSZOwdetWnDt3Di9evIC/vz+CgoIqnCQVFBSQmZmJtLQ0lJSUoLCwEBkZGXy9FL81fvx4bNu2DRcvXkRqaiqWLFkCJSUldOjQ4af7atKkCeLi4pCcnIykpCR4eXkhOTmZrylS2Dw9PZGcnIxp06YhPj4eL168QEhICGbMmAFLS0v07dsXcXFx2L17N6ZMmYJGjRohIyOD7/W9XmanT5+GmpoaRowYwfd71aNHD5iamiI0NBSlpaVYtWoVSkpK4O3tjaKiIl69Hz9+FNnnUN1R85qEOXPmDAYPHsx3H+IrR0dHHDp0CLKysmjQoAG6dOkCe3t7KCgoYNKkSZgwYQIAQE1NDTt37sT69etha2uL+vXrw8bGhu8Bux9p2rQpdu7cCR8fH9ja2kJFRQV2dnZwd3cHADRv3hw7d+7E6tWrceDAAejp6cHe3h6PHz8uU5e0tDS2bduGlStXYtSoUahXrx6GDx+O6dOnl7tvHx8frFy5EuPGjYOcnBzMzc2xYMECAICHhwc4HA4WLlyIvLw8dOnSBfv37+c1i1TFkCFDEBERgXXr1qF9+/bYvXs3736Dk5MT8vPzsW7dOrx//x7a2trYvn17mftA3+Pi4oLCwkL4+PggMzMTrVu3hp+fX7ldf8tjYWGB48ePY/DgwTh06BCSk5OxcOFCREZG8t3P+2rixInIy8uDl5cX8vLyYGBggICAgB82pX61ePFiLFmyBCNGjEC9evXQs2dPTJ06VejPh31LU1MTR48exdatW+Hh4YGcnBxoaGhgwoQJGD9+PDgcDm7cuAHGGPz9/eHv71+mjuHDh8Pb25uv7Oszby4uLuXem3FycoKrqyuio6MRGRkJLpeLIUOG8K1jYWHB192+NuOwql7/khorODgYa9euRVxcnLhDqdGcnZ3RuXNnXlKrCaZMmYLVq1d/txMIIcJGzWuE1BK3b9+GlJQUJRwiVpR0CKkldHV1sW3bNnGHQWo5al4jhBAiMnSlQwghRGRqddIpKSlBenp6lZ7gJoQQUnm1Oum8efMGAwcOxJs3b8QdCiGE1Aq1OukQQggRLUo6hBBCRIaSDiGEEJGhpEMIIURkKOkQQggRGUo6hBBCRIaSDiGEEJGhpEMIIURkKOkA353UihBCiGBR0gHw5MkTcYdACCG1AiUdfJlnhBBCiPBR0gElHUIIERVKOgDu3r1LI00TQogIUNIBkJeXh4cPH4o7DEIIkXiUdP4nNjZW3CEQQojEo6TzP9euXRN3CIQQIvEo6fzPjRs36L4OIYQIGSUdAM2bN8enT5/ovg4hhAgZJR0ABgYGAOi+DiGECBslHfx/0qH7OoQQIlyUdMB/pVNQUCDmaAghRHJR0gHQpEkTdOnSBfn5+bh8+bK4wyGEEIlFSed/LC0tAQAXLlwQcySEECK5KOn8j5WVFQAgPDycuk4TQoiQUNL5H21tbWhqaiIrKwvx8fHiDocQQiQSJZ3/4XA4vKud8+fPizkaQgiRTJR0vmFhYQHgy30dxpiYoyGEEMkjUUnnwIEDCAwMrPL2hoaGUFVVRVpaGo1OQAghQiAxSWf58uXYv3//L9UhJSUFc3NzANSLjRBChEFiko6JiQmmTZv2y/V8va9z9uxZamIjhBABk5ikY2ZmJpB6TExM0LBhQzx58oSa2AghRMAkJukIiqysLIYOHQoACAkJEXM0hBAiWSjplGP48OEAgNDQUJSWloo5GkIIkRyUdMrRrVs3tGzZEm/evKHpDgghRIAkKunY2dlhypQpv1wPh8PhXe1QExshhAhOtU069+/fh5GREV/Z48ePYW9vDz09PVhaWiImJkZo+7ezswMAhIWFobCwUGj7IYSQ2qRaJp1z585h4sSJKC4u5pUVFRVh2rRpsLCwwM2bNzF//nzMmTMHL1++FEoMbdu2RZcuXfDp0ydEREQIZR+EEFLbVLuk4+vri507d2L69Ol85XFxcSgsLMT48eMhKyuLgQMHwsjICGfOnBFaLNTERgghglXtko6DgwOCg4Oho6PDV/706VNoaWmBw+Hwytq0aYOkpCShxWJrawsOh4OoqCjk5OQIbT+EEFJbVLuk07Rp03LL8/PzIS8vz1cmLy8v1Oml1dTU0KtXLxQVFdHI04QQIgDVLul8j4KCQpkb+oWFhVBQUBDqfr92KKAmNkII+XU1JuloaWkhNTWVrywlJQVt27YV6n6trKwgJyeHa9eu4c2bN0LdFyGESLoak3SMjY0hLS2NwMBAFBUVISoqCnFxcRgyZIhQ91u/fn0MHDgQjDGcPn1aqPsihBBJV2OSjpycHHbs2IGLFy+iR48e8PHxga+vL1q0aCH0fQ8bNgzAl2FxCCGEVB2H1eLx+9PT0zFw4EBERkaiefPm312voKAAenp6yM3NxaVLl6ClpSXCKAkhRHLUmCsdcapbty5vnh3qUEAIIVVHSaeCRowYAQA4ceIEuFyumKMhhJCaiZJOBfXu3Rvq6upIS0tDXFycuMMhhJAaiZJOBUlJSWHkyJEAgGPHjok5GkIIqZko6VTCqFGjAABnz55FXl6emKMhhJCah5JOJbRp0wbdunVDfn4+zp07J+5wCCGkxqGkU0mjR48GQE1shBBSFZR0Kmno0KGQl5fHtWvXkJaWJu5wCCGkRqGkU0nKysqwtLQE8KX7NCGEkIqrctLhcrlITk5GXFwc4uPj8ezZM9SWwQ2+dig4ceJErTlmQggRBJnKbhAfH4/9+/fj6tWrfD24OBwOlJWV0adPHzg4OMDAwECggVYnJiYmaNq0KZ49e4bbt2/D0NBQ3CERQkiNUOGk8/z5cyxduhSvXr2CmZkZ/P39oaWlhQYNGoDL5SI7Oxv//vsvbt68CQ8PD7Ro0QJ//PEHWrduLcz4xUJGRgbDhw/H9u3bceLECUo6hBBSQRUe8HP06NFwd3eHqanpT9dljCEyMpL3pVxdVXTAz/I8evQIgwYNQoMGDXD79m3UqVNHSFESQojkqPA9naCgoAolHOBLU5uZmRmOHz9e5cCqu06dOqFTp07IyclBZGSkuMMhhJAaocJJh8Ph8P6/aNEiPH78uMw6WVlZGDhwYLnbSKKvw+JU56s5QgipTqrUey04OBj29vY4e/YsXzmXy8WrV68EElhNMHz4cEhJSSEyMhLv378XdziEEFLtVbnL9Jw5c7BkyRKsWbOm1g7136RJE/Tr1w8lJSXYsGGDuMMhhJBqr8pJx9raGocPH8Y///yD8ePHIysrS+Kb08qzcOFCyMjIYN++fbh165a4wyGEkGqtSknna3Lp2LEj736GnZ0dnjx5IrjIaohOnTph6tSpYIxhwYIFKC4uFndIhBBSbVUp6Xzby7phw4bYs2cPzMzMMHXqVIEFVpPMnTsXLVu2RGJiIgIDA8UdDiGEVFtVSjozZsyAgoIC7720tDS8vLywYsUKdOvWTWDB1RR169aFt7c3AGDDhg1ISUkRc0SEEFI9Vfjh0OTk5ApX2rZt2yoHJEq/8nBoeWbOnIng4GAYGhoiODgYMjKVHmWIEEIkWoW/Fa2trcHhcMAY4/sX4G9u43A4SExMFHykNcDKlStx7do13Lp1C9u2bcPMmTPFHRIhhFQrFU46/33qnjGGoUOHIjAwEOrq6gIPrCZq0KABNmzYAAcHB/z1118YMGAAdHR0xB0WIYRUGxVOOhoaGuWWq6mpfXdZbWRqagoXFxfs378fs2fPxrlz5yAnJyfusAghpFqgSdyEYMmSJdDU1ERiYiK2bt0q7nB4rl+/Dh0dHezfv1/coRBCailKOkKgoKCAdevWAQA2bdpUqU4YwlJaWorFixcjJycHq1atwrt378QdEiGkFqKkIyS9e/fG2LFjUVRUhPnz54t9qKBjx47h33//BQDk5eXBx8dHrPEQQmqnCt/TOXToUJmy0tJSnDp1CioqKnzljo6Ovx6ZBPDy8kJERARu3LiBQ4cOwdnZWSxx5Ofn85LMvHnzsGnTJhw5cgTjx4+njg6EEJGq8HM6AwYMqFiFHE6NmV9G0M/plOfMmTNwc3ODsrIyrl+/jvr16wtlPz/i6+uL9evXQ09PD2fPnsWyZcuwa9cumJiY4OjRo7VyzDxCiHhUOOlIIlEkHcYYRo0ahdjYWMyfPx9z5swRyn6+Jzk5GVZWVsjPz8fJkyfRo0cPZGdnw8TEBDk5Odi5cyesrKxEGhMhpPaq8D2dQYMGYcWKFYiIiEBubq4wY5IoHA4HHh4eAIAdO3aI9LO7efMmbG1tkZ+fj8GDB6NHjx4AABUVFcybNw8A4OnpiczMTJHFRAip3SqcdAIDA9G2bVucPHkSAwYMgL29Pfz9/XH79m2x3ySv7nr27AkjIyPk5ORg7969ZZYnJiZi0qRJiIqKEtg+z5w5gzFjxiAnJwcDBw7Exo0b+ZaPHz8evXr1QmZmJubPn49afMFLCBGhKjWvlZaW4s6dO7h27RquXr2KZ8+eoVu3bujduzd69+6NVq1aCSNWgRNF89pXly5dwtixY9GwYUNcv34dioqKAID4+Hi4uLjgw4cPkJeXx4kTJ6Cvr1/l/cTHx+Ovv/7CpUuXAABOTk74888/yx0H7uXLlzAzM8PHjx/h4+MDBweHKu+XEEIqQiD3dD59+oTY2FhcvXoVsbGxCA8PF0RsQifKpMMYg42NDW7fvg1nZ2dYW1vj48ePmD17NvLz86GmpoY3b95AVVUVYWFh0NDQwOPHj/Ho0SNYWVlBXl6eV1deXh6Sk5OhoaGBRo0aISsrC+fPn0doaChiY2MBAEpKSpg3bx4mT578w44CISEhvFHDg4KCYGBgINTPgRBSu1FHAhElHQCIiooqt9v0yJEj4e3tjXHjxuHq1avQ0tJC3bp18eDBAwBA9+7dsXv3bjRs2BA3b97ElClTeA93KikpoaCgAKWlpQAARUVFuLq6YsqUKWW6sn/PjBkzEBISAllZWSxZsgQTJ07kS1SlpaU4ePAgXrx4AWNjY/To0QPKysq/+nEQQmqhCiWdp0+f4t27d9DV1eU1CwFAdHQ0+vfvL9QAhUnUSYcxhgMHDiA+Ph4vXrzAu3fvYGNjg99//x1SUlLIycnB0KFDefPxKCsro06dOsjIyICmpiZGjRoFX19flJSUoFmzZsjNzcWnT58gIyODPn36wNraGpaWlmjQoEGl4ioqKsKqVauwa9cuAIClpSXmzJmDLl264N27d3B3d8e1a9d460tJSaFFixZo1KgRGjduDGdn5wp3qSeE1G4/TTr79+/H4cOH0aZNGyQmJmLhwoUwNzcHAAwfPhwhISEiCVQYRJ10KiItLQ07d+5E9+7dYWZmhpycHIwfPx4JCQm8dSZPngwvLy9IS0sjOzsbsrKyqFev3i/v++zZs5g3bx6vh13nzp3x9u1bZGRkoHHjxhg5ciRu3bqFO3fuoKSkhLddq1at+JISIYR8F/sJa2trlpeXxxhjLC0tjQ0fPpzt3r2bMcaYra3tzzav1tLS0pi2tjZLS0sTdyg/lJeXxyZPnsw6dOjATp48KdR9PX/+nC1ZsoR16tSJqaurM3V1dTZy5Ej25s0b3jr5+fns6dOnLC4ujmlrazN1dXW+5YQQ8j0/vdKxtrbG2bNnee/z8vIwa9YstG3bFtevX8epU6eEnhiFpTpe6fwIl8uFlJRohssrLCzEhQsXkJeXhzFjxkBaWrrc9caOHYtLly4hICAA1tbWIomNEFJz/fQbrGHDhnwzgSoqKiIgIADZ2dl48uSJUIMj/ESVcABAXl4etra2cHBw+G7CAQAjIyMAwI0bN0QVGiGkBvvpgJ/r1q0r86UjIyODdevWYcyYMUILjNQM3bp1A/Dl+SBCCPmZn/7prKamBlVVVQAo05RmaGgonKhIjWFgYABpaWk8ePAAeXl54g6HEFLNVaq9ZtGiRXBxccHTp0+FFQ+pYRQVFaGjo8MbpYIQQn6kUknn5MmTKC0tha2tLXx8fFBQUCCsuEgN0r17dwBfBhglhJAfqVTS6dChAw4dOoQ///wTp0+fhpWVFS5cuCCs2EgN8TXpUGcCQsjPVKk7lK2tLf7++29YWVlh3rx5mDRpEl68eCHo2EgN8TXp3Lp1i++hUUII+a8q98FVVFTEggULcOrUKXC5XFhbW8PPzw9FRUWCjI/UAGpqamjZsiXy8vLw77//ijscQkg19ssPfmhpaWH37t3466+/EBoaiiFDhggiLlLDUBMbIaQifvqczve8fPkST548wZMnT/D48WM8efIEb9++pQndaqnu3bvj5MmTiIiIwMSJE8UdDiGkmqpU0lm6dCmePHmCpKQk5OfnAwBatGiBDh06wMrKCh06dECHDh2EEiip3qysrLBy5UrExMTg8uXL6NOnj7hDIoRUQ5VKOo8fP0b79u1hY2ODDh06oH379nxTHZDaq3HjxpgxYwbWrl2LP/74A3///fcPh88hhNRONIlbDRrws7orKChA37598erVK/z1118YM2YMzp8/j6ioKMyZM4c+Y0JI1ToShIaGIikpie7fED5169bFwoULAQBr167F8OHDMXnyZBw5cgRz5sxBLf77hhDyP1W60unRowc+fPgAOTk5tGvXDp06dULHjh3RqVMntG/fHvLy8sKIVeDoSkfwuFwuhg4dirt37wIAGjVqhNLSUuTk5MDPzw8jRowQb4CEELGqcvPa69ev8fDhQyQmJuLevXu4fv06SktLISUlhYcPHwo6TqGgpCMcDx48gKenJ/r06YPp06fj3Llz8PDwQOPGjXHp0iXUr19f3CESQsREYPd03rx5g7lz58LW1hb29vaCqFLoKOmIBpfLxYgRI3Djxg2MHz8ef/75p7hDErhTp04hNjYWDg4O0NXVFXc4hFRbApsVTE1NDWvWrEFMTIygqiQSQkpKCmvWrIG0tDT27t2Lli1bQkNDA8bGxrxmOAD4+PEjJkyYAGdnZ2RkZIgv4Ep6/PgxZs+ejQMHDsDKygpjxoyh+YUI+Y4qJZ2srKxyyzU1NfHo0aNfCohIpg4dOmDOnDkAgNLSUgBfrjQdHBzw4MEDfPz4EY6OjggPD0dUVBSGDh1aI2amLS0txbx581BcXAxdXV0oKiriypUrsLe3/+7vCSG1WZWSTq9evWBqago3Nzf4+fkhIiICKSkpCA0NRXFxsaBjJBLCw8MDqampePbsGVJSUmBpaYkPHz7A3t4eY8aMwe3bt9GiRQt07doVaWlpsLW1xZUrV8QdNj58+ICAgIByk8ju3btx584dqKmpISgoCDdv3oS+vj4KCgpw+fJlMURLSPVWpaRz5MgRTJ06Faqqqrh48SI8PDwwePBgLF68GFOnThV0jESCyMnJQVZWFnXq1MHWrVsxYMAAZGdn4/79+2jRogWOHz+OkydPYvDgwfj48SMmTZqEt2/fii1exhjmzp2LP/74A+7u7nzdvl+8eIG1a9cCANasWQNlZWXUr18fgwcPBgBcunRJLDETUq0xASgpKWHPnz9nOTk5gqhOZNLS0pi2tjZLS0sTdyi1VkFBAXNzc2ODBw9mL1684JWXlpYyJycnpq6uzqZNmya2+EJCQpi6ujrvdfToUcYYY3l5eWzIkCFMXV2dubm58W2TkJDA1NXVWbdu3RiXyxVH2IRUWxXuvfbx40coKytXKqFVZRtRot5r1VtaWhr69euHwsJCHD58GKampgKre8eOHZCRkcGECRO+u05mZib69euH7OxsmJubIzw8HPXr18c///yDBQsWIDo6Gs2bN0dYWBgaN27M247L5UJPTw/v379HTEwM2rZtK7C4CanpKty85uTkhMDAQHz69Omn62ZlZWHr1q1wdHT8peBI7daiRQt4eHgAABYuXCiw6dHj4+OxfPlyeHl5ITw8/LvrLV68GNnZ2ejbty92796N/v3748OHDzA3N0d0dDQaNmyIw4cP8yUc4Etvva8DnlITGyH8Knylk5eXh40bN+LEiRMwMDBAnz590K5dO6ioqIDL5SI7Oxv//vsvbt68iZs3b8LOzg5z586FkpKSsI+hyuhKp/orLi6GpaUl/v33X+jp6aFjx45o3rw5nJ2d0ahRoyrV6ezsjKioKABA06ZNERUVhQYNGgAAnj59itOnTyMsLAyJiYlQVFREVFQUmjdvjvT0dAwYMAB5eXlQUFDA8ePHoaenV+4+goKC4OHhgUGDBmHv3r1VipMQSVTph0MzMjJw7NgxXLx4EY8ePeJ1f5WRkUGnTp1gamqKkSNHomnTpkIJWJAo6dQMN2/ehL29PQoLC3llXbp0QVhYWKVHsr579y6GDBkCBQUFaGtr4+7duxg5ciT+/PNPeHt7Y+/evbzOAsrKyli7di1sbGx42586dQr+/v5YsmQJ+vbt+939vHr1Ct27d4eioiIePnwIWVnZSh41IZLpl0Yk4HK5yMnJAYfDgYqKiiDjEglKOjXHq1ev8OjRI7x69Qp+fn54/fo11q5dCycnp0rVM2HCBISHh2PatGkYO3YszM3NUVhYiMaNGyMzMxMyMjIYNmwYbGxs0KdPH8jJyVU55n79+iEpKQnBwcEwNjZGfn4+5OTkICNT5bkTCanxfmlEAikpKTRs2LBGJhxSs6irq8PMzAwuLi5YtmwZgC/dlCvzAOaDBw8QHh4OeXl5TJ06FVpaWliwYAGAL50G9PT0cP78eWzatAkDBw78pYQDgHcldOLECSxcuBA6OjqwtrZGZmbmL9VLSE1WpT+5OnToAA6HU+4yWVlZNG3aFIMHD8asWbNoIi8icNbW1ujVqxeuXbsGHx8frFmz5qfbFBQUYPXq1QAAR0dHqKqqAgBcXV1RWFgIFRUVODg4CPTntW/fvti1axcOHz7MK0tISMCIESNw9OhRNGvWTGD7IqSmqFLz2rFjx7Bp0ybMnDmTdyM1ISEBfn5+GD16NNq1awd/f38MHDgQc+fOFXTMAkPNazXX48ePMWjQIHC5XJiZmSErKwufP3+Gubk5xo4dC3V1dd66Dx8+hLu7O5KSkqCoqIiYmBiRfOHn5eWhZ8+eyMnJwbBhw+Do6IjFixcjMTERLVu2xNGjR9GqVasy2xUUFCA8PBxJSUlwc3Pj64xTUlKC0tJS1KlTR+jxEyIUVXm4x9zcnEVHR5cpv3TpErOwsGCMMRYfH8/69u1blepFhh4OrdmWLVvG9+Dm11fz5s3ZiBEjmIuLC3NycmKamppMXV2d9e3blyUkJIg0xtevX7PXr1/z3mdlZbHBgwczdXV1pqOjw/d7lJCQwGbPns3atWvHO5YlS5bwlpeWlrIxY8awNm3asD179pT74CmXy2VhYWHs3r17ZcqfPXtGD6sSsatS0tHT02PJycllyp88ecJ0dXUZY4ylp6ezrl27/lJwwkZJp2b7/PkzCwkJYWFhYSwuLo5FR0ezqVOnslatWpVJRL///jvLz88Xd8iMMcY+fvzIG21BQ0OD/fHHH2z8+PF88VpaWjJ1dXWmqanJ0tPTGWOMBQcH863j5OTE3r17x6u3tLSULV68mLd83rx5LDMzk126dIlZW1szdXV15unpKa7DJoQxVokRCb7l6uoKWVlZeHt7855vyMnJgaenJz5//ow9e/Zg7969CA0NRWhoqICvzQSHmtckU2ZmJm7fvg0A4HA4aNasGTp37izmqPhxuVxs2rQJf/31F6+Ltry8PJydneHi4oI2bdpg2rRpOH36NBwdHbFixQr07dsXr169wtixY3H+/Hnk5OSgfv36mDRpEiZMmIDVq1fj8OHDkJOTA2MMxcXFqFOnDj5//sy370OHDqFfv35iOGpCqnhP59WrV5gyZQpevHgBDQ0NcLlcvH79Gm3atIG/vz9SUlIwY8YM+Pn5Vesfbko6RNxiYmLg4+ODHj16wM3NjW90g+TkZPTv3x8cDgejR4/GkSNHoKOjg/Pnz+Pdu3eYP38+oqOjAXzpwFNcXAx5eXns3r0bzZs3h5eXF2+mVjc3NxQVFcHX1xfNmjVDVFRUtR6iikiuKj+nU1paiuvXr+Px48eQkZFBu3bt0LNnTwD/P99Ow4YNBRepEFDSIdWdh4cHgoKCeO+PHTuG3r17897HxsZiw4YNuHbtGhQUFLB//37e7yFjDAkJCWjVqhXq16+PkpISDBs2DHfu3IG9vT3++usvkR8PIVV+TkdaWhpSUlKQkZEBl8tFcXExSkpKAHxJNtU94RBSE8ydO5c3moGFhQVfwgGAnj174vjx4/j7778RERHBSzjAl6ZFXV1d1K9fH8CXUUN8fX1Rp04dHD16lMaFI2JRped0MjIy4O7ujkePHvE1r7Vu3Rp79+6t8phYhBB+LVq0wNy5c3Hy5EksXbr0u+t16dKlQvW1a9cOs2fPxrp16xAQEPDDoXwIEYYqNa/NmTMHb9++xcaNG3ljrL158wbz5s1Ds2bNsH79eoEHKgzUvEZqo+zsbBgaGqKoqAjXrl1Dy5YtxR0SqUWq1Lx2+fJleHl58Q3qqaamBk9PT7pkJ6SaU1FRwZAhQ8AY4xstISsrCzdv3hRjZKQ2qFLSqVOnDqSkym4qJSXFG3WaEFJ9fR0oNSgoCMXFxfjw4QNsbGwwbNgw+sORCFWVkk6fPn2wevVqvoELMzIy4O3tDRMTE4EFRwgRDiMjI7Rr1w7v3r3DhQsXMHPmTKSmpgIAtmzZIuboiCSrUtL5/fffkZ2djf79+8PCwgIWFhYYMGAACgoKsHjxYkHHSAgRMA6Hw7va+e233xAZGYkGDRpAUVERV65cwYMHD8QcIZFUv/SczqVLl/D06VPUqVMHWlpa6NWrl6DjEyrqSEBqs68dCj5//gwpKSkcPHgQ0dHR2LFjB+zs7LB582Zxh0gkUIWTTnJycoUrbdu2bZUDEiVKOqS2+/3333Ho0CEsWrQI7u7uSE9P5/3xGBsbCw0NDTFHSCRNhZPO1zl0vl39v3PqMMbA4XCQmJgo2CiFhJIOqe0KCwvx9OlT6Ojo8Mrc3d0RGhqKqVOn/vDZIEKqosJJ5+XLlxWuVJR/HZWUlGDBggV4/fo1VFVVsW7dugrPNUJJh5Cy7t+/DysrKygpKeHvv/9G69atxR0SkSAV7kigoaFR4ZcohYeHo2nTpjh8+DB0dXUREhIi0v0TIml0dXVhbm6O3NxcjB49GmlpaeIOiUiQKo+9Vl3cvXuX1wbdu3dv3LhxQ8wREVLz+fv7o3v37nj16hXGjBmDV69eiTskIiFqfNLJzc2FoqIiAEBRURF5eXlijoiQmk9RURH79+9H165d8fz5c4waNarG3Ksl1VuNTzpKSkrIz88H8GVO+m/nkyeEVJ2ysjIOHTqELl264NmzZ7C2tuabZoGQqqj0KNNZWVm4ePEikpKSeF/yHTt2hKmpqVgmhdLV1cW1a9fQu3dvXL16Ffr6+iKPgRBJpaKigpCQECxevBhBQUHw8PDAgwcP8Mcff5TpvUpIRVTq4dCdO3fC398fHA4HGhoaUFZWRm5uLlJTUyEjI4O5c+fCxcVFmPGWUVJSAk9PT6Snp0NFRQW+vr6Ql5ev0LbUe42QigsKCsKiRYtQWFiItWvX8kY0IKQyKnylc/z4cfj7+2P+/Pmws7ND3bp1ecsKCwsRHByM9evXo2nTprCwsBBKsOWRkZGpMVMpEFKTjRkzBrKyspg5cyaWLl2Krl27VngeH0K+qvA9nYMHD+K3336Do6MjX8IBAHl5eTg4OGDu3LnYv3//Lwd1//59GBkZ8ZU9fvwY9vb20NPTg6WlJWJiYn55P4SQyrGzs4OTkxM+f/6MqVOnIjY2Ftu3b4e7uzsWLlyIHTt2ICYmBlwuV9yhkmqqwlc6z58//+kI0iYmJr88Qu25c+ewdOlSvikSioqKMG3aNDg7O+PAgQO4dOkS5syZg7Nnz9IwHYSI2IoVK3Dv3j0kJCRg5MiR5a5jbW2N7du3030fUkaFr3QKCwt/2jOsXr16+PDhQ5WD8fX1xc6dOzF9+nS+8ri4OBQWFmL8+PGQlZXFwIEDYWRkhDNnzlR5X4SQqpGXl0dAQABat26Ntm3bwsHBAT4+Pli+fDnGjRuHevXq4ezZs9i4caO4QyXVUKV6rwn7r5avTXRxcXF85U+fPoWWlhbf/tu0aYOkpCShxkMIKV+rVq1w5cqVcpcNHDgQ48aNw/r169GhQwdYWVmJODpSnVUq6QQHB0NBQeG7y3/1wcxvp7/+Vn5+fpkeafLy8igoKPil/RFCBG/gwIFYvHgxVq1ahZkzZ6J3795QVlaGuro6Zs+e/cPvECL5Kpx01NXVceTIkZ+u16xZs18KqDwKCgooLCzkKyssLKQfXkKqKTc3N/z77784ceIEIiIieOUcDgeenp5ijIyIW4WTTlRUlDDj+CEtLS3s3LmTrywlJYUeBCWkmuJwONi4cSNcXFyQmZmJFy9eYPny5dixYwcmTJjw3VYNIvlqxDA4xsbGkJaWRmBgIIqKihAVFYW4uDgMGTJE3KERQr6Dw+HA0NAQFhYWmDx5MgYPHozCwkJs2LBB3KERMapU0snNzUVoaChyc3MBfJm0bceOHXBxccGcOXNw//59oQQpJyeHHTt24OLFi+jRowd8fHzg6+uLFi1aCGV/hBDBW7BgAaSkpHDkyBEkJyejqKgIoaGh2Lt3LzIyMsQdHhGRCg+D8+LFCzg6OiIvLw+nTp1CixYtsGrVKhw6dIg34dOZM2ewa9cuGBgYCDtugaBhcAgRra/TY3fu3BmZmZl48+YNgC8ji1hYWMDZ2RkmJib0fI8Eq3DSmTdvHvLy8rBx40bIy8vj/fv36Nu3L8zNzeHr6wsA2LVrF65evYrdu3cLNWhBoaRDiGi9fv0aJiYmvI5B2traaNmyJaKionijGHTs2BGTJk3CsGHDKjyOIqk5Kty8FhsbC3d3d94PweXLl8HlcjFs2DDeOr1798a9e/cEHiQhRDI0a9YMa9euxeDBg7Fv3z5ERkZi3759uHHjBubNmwdVVVUkJiZi3rx5cHR0RCXGIyY1RIWTzqdPn9CoUSPe+7i4OMjIyPCNkaaoqEhjLhFCfmjkyJHYsWMHzMzMICX15SuoWbNm8PDwQFxcHDZt2oSGDRvi+vXrfN2tiWSocNLR0NBASkoKAKC0tBSXL19G9+7d+Qb/vH79OjVTEUKqrE6dOhg5ciRmzpwJANi4cSNd7UiYCicdOzs7rFy5EmfOnMGSJUuQmZmJsWPH8pbHx8dj48aN1I2ZEPLLnJ2d0ahRI9y9e5dGlJcwFU46rq6uGDhwIFatWoXo6Gh4eHhg0KBBAIBVq1bByckJhoaGcHV1FVqwhJDaoW7dupg6dSoAutqRNJWaOfR7Hj9+DC6Xi44dOwoiJpGh3muEVF+5ubkwNjZGTk4ONm/ejK5du0JJSQlNmjShLtU1WIWHwcnPz8eaNWtw4cIFyMnJwczMDL/99huUlJTQvn17YcZICKmFlJSUMHnyZPj4+PDu8QCAmpoa+vfvj4EDB2LQoEGQkanUuMVEzCp8tjZv3ozo6GhMmjQJ0tLSOHToEHJycmjODEKI0EyaNAkJCQl4/vw5CgoKkJOTgzdv3uDIkSM4cuQIzM3NERgYCFlZWXGHSiqowknn/Pnz8PHxQc+ePQEA3bt3h6OjI4qLi+mEE0KEQklJCbt27eK9Z4zh4cOHiI6Oxvbt2xEeHo4ZM2Zgy5YtdMVTQ1S4I0FGRga0tLR477t06QLGGN6/fy+UwAgh5L84HA46d+6MmTNn4vDhw7xZSufOncs3xT2pviqcdEpLSyEtLc17z+FwICsri+LiYqEERgghP9K1a1ccPHgQCgoKCA4Ohr29PV6+fCnusMhP0PUoIaTG6tatGw4ePIgpU6bg2rVrGDRoEObOnYvMzEzcunUL2dnZ6NSpE/T09NC/f3+0adNG3CHXer80XXVpaSlOnToFFRUVvvUcHR0FEx0hhPyEsbExIiMj4eHhgcjISCxfvpxv+b///ovg4GDIycnh+PHj6Natm3gCJQAq8ZzOgAEDKlYhh4PIyMhfCkpU6DkdQiQHYwyHDx9GZGQk2rVrB0NDQzRu3BgJCQn4559/EB0dDTU1NZw/fx5NmjQRd7i1lkAeDq2pKOkQUjsUFxfD3t4e169fh5GREYKCgiAnJyfusGqlGjFdNSGE/ApZWVls374dampquHHjBjw9PfH582dxh1UrUdIhhNQKqqqq2LFjB+Tk5BAUFARLS0vcunVL3GHVOpR0CCG1hoGBAY4dO4bWrVvjyZMnsLW1xfz58/H48WNxh1ZrUNIhhNQq3bt3xz///IMZM2ZASkoKhw8fxoABAzB69GhcvXpV3OFJPEo6hJBap27duli4cCGioqIwbtw4KCgo4OrVqxg9ejSmT5+O169fiztEiUVJhxBSa7Vt2xarV69GfHw8fvvtN8jLy+PUqVMwNTVFWFiYuMOTSJR0CCG1Xv369TF37lzExMTA0tISeXl5mDVrFhISEsQdmsShpEMIIf/TvHlz7Ny5Ew4ODigsLISrqysNaixglHQIIeQbHA4Hq1atgoGBAV6+fImpU6fSMz0CREmHEEL+o06dOtixYweaNGmC2NhYdO7cGS4uLtizZw8+ffok7vBqNEo6hBBSDjU1Nezduxc6OjrIz89HZGQkvLy80KtXL+zYsQOFhYXiDrFGorHXaOw1QshPvHnzBpcuXcKhQ4cQHx8PAGjWrBmcnJxgb28PNTU1MUdYc1DSoaRDCKkgxhgiIiKwdu1aJCYmAgCkpaVhaWmJVatW0ejVFUDNa4QQUkEcDgeDBg1CeHg4jhw5giFDhoDD4SAsLAwWFha4du2auEOs9ijpEEJIJUlJSaFv374IDAxEbGwsevbsiXfv3mHMmDFYuXIlbt68ieLiYnGHWS1R0iGEkF+grq6Oo0ePYtasWeByudi+fTuGDRuGTp06YeLEiThz5gwKCgrEHWa1Qfd06J4OIURA4uLiEBoaiqtXr+Lp06e8ciUlJYwYMQLTpk1DixYtxBih+FHSoaRDCBGC169f4+zZswgJCcG9e/cAADIyMhg+fDhmzZqFNm3aiDlC8aDmNUIIEYJmzZph8uTJOHfuHCIjI2FnZwfGGI4fP45+/fphwYIFePPmjbjDFDlKOoQQImQdOnTA5s2bcfnyZYwdOxaMMRw8eBC9e/fGnj17xB2eSFHSIYQQEWnVqhXWr1+P6OhoDBkyBIWFhfDy8oKXlxdKSkrEHZ5IUNIhhBARa9u2LQIDA7F582bIyclhz549mDBhAhITEyHpt9kp6RBCiJjY2dkhKCgIKioqiIqKgpmZGQwNDTF//nxkZWWJOzyhoKRDCCFiZGRkhLCwMIwePRpNmzbF27dvcfjwYTg7OyMvL0/c4QkcJR1CCBGzVq1awdfXF7du3cKFCxfQsmVL3L17F5MnT0ZRUZG4wxMoSjqEEFJNcDgcdO7cGYcOHUKjRo0QExOD2bNnIzk5WWLu9VDSIYSQaqZNmzY4cOAAFBQUcPr0aZiamsLQ0BAeHh5ITk4Wd3i/hJIOIYRUQ127dsWxY8dgbW2Nxo0b4+3btwgKCkL//v0xd+5cpKWliTvEKqGkQwgh1ZS+vj4CAgJw9+5dREZGwtHRERwOB8eOHYOFhQWeP39eZhsul4vw8HDMnDkTFy9eFH3QP0Fjr9HYa4SQGiQ1NRULFizA1atX0bVrV4SGhkJOTg6lpaU4ePAgAgMD8ezZMwBA/fr1ERMTA1VVVfEG/Q260iGEkBqkdevWCAwMRPPmzXHv3j2sWbMGGRkZcHR0xKJFi/Ds2TM0b94cOjo6+PDhA5YvXy7ukPlQ0iGEkBqmQYMG2LJlC6SlpREYGIj+/fvj8uXLaNSoEbZu3YqrV69i586dkJeXR2hoKKKjo8UdMg8lHUIIqYG6deuGBQsWAACys7PRs2dPhIeHw9bWFjIyMmjZsiV+++03AMCiRYvw5s0bfP78+Yd1FhQUYMOGDVi5ciVKS0uFEreMUGolhBAidNOmTUNRUREUFBTg6uoKGRn+r/RJkyYhODgYjx49gqGhIQCgXr16GD16NNzd3dG0aVPeujdv3oSHhwdSUlIAAMbGxjA3Nxd4zNSRgDoSEEIkWGJiIubNm4fnz5/j06dPvCsYeXl52NjYoKioCOnp6bh16xYYY1BSUkJubi7MzMywb98+gcdDSYeSDiGklmCM4dGjR/D19cX58+f5lklLS2P69OlwcXFBr169UFpaiuvXr0NDQ0OgMVDzGiGE1BIcDgc6OjrYuXMnHjx4gMuXL6Nx48bQ0NBA27Zt0aRJEwCApaUlzpw5g6CgIHh4eAg2BrrSoSsdQgj51uXLl2Fvbw91dXVcv34d0tLSAqubeq8RQgjh07t3b7Rq1QqvXr0S+KgGlHQIIYTwkZKSgoODAwDgwIEDfCNcM8Zw8+bNqtf9y9ERQgiROKNHj4asrCz++ecfjB8/Hq9fv0ZycjLGjBmDYcOGVble6khACCGkjCZNmmDTpk3w9PREREQE+vXrh8+fP6O4uBgqKipVrpeudAghhJTL1tYWUVFRMDMzQ25uLoqLi+Hg4IBLly5VuU660iGEEPJdzZo1w969exETE4OGDRtCV1f3l+qjpEMIIeSHOBwO+vXrJ5C6qHmNEEKIyFDSIYQQIjKUdAghhIgMJR1CCCEiQ0mHEEKIyFDSIYQQIjK1usv018mM3rx5I+ZICCGk5lFTUyszW+nP1Oqkk5GRAQBwdHQUcySEEFLzVGVamFo9n05hYSEePHgAVVVVgc4XQQghtUFVrnRqddIhhBAiWtSRgBBCiMhQ0iGEECIylHQIIYSIDCUdQgghIkNJhxBCiMhQ0iGEECIylHQIIYSIDCWdGubYsWPQ0dGBvr4+7xUSEiLusATm/v37MDIy4r0vKirCsmXLYGxsDGNjY6xZs4Y3fFFN9t/jLCwsLHNeJ06cKMYIq+7+/ftwcHBAt27d0K9fP2zevBmMMYk7l987Tkk6lwAQHR2NoUOHQl9fH4MGDcLRo0cB/MLvJiM1ypIlS9iGDRvEHYZQhIWFMUNDQ6anp8crW79+PXNwcGDZ2dns9evXzNbWlgUEBIgxyl9X3nHeuXOH9enTR4xRCUZubi7r0aMHO3DgACspKWGpqalswIAB7ODBgxJ1Ln90nJJyLhlj7OXLl6xz587s8uXLjDHGEhMTWZcuXdi9e/eqfD7pSqeGefjwITp27CjuMATO19cXO3fuxPTp0/nKQ0JCMGXKFDRo0ABqampwc3PDyZMnxRTlr/vecUrKeX39+jUMDAzg5OQEaWlpaGpqYtCgQbh9+7ZEncsfHaeknEsAUFdXR2xsLExMTMDlcpGdnQ1paWkoKSlV+XxS0qlBiouL8eTJEwQHB8PExASDBg1CYGAgmASMZOTg4IDg4GDo6Ojwyj5+/IiMjAxoaWnxytq0aYPnz5+jqKhIHGH+svKOE/iSdDIyMjB06FD06tULs2bNwtu3b8UUZdW1bdsWW7Zs4b0vKirCpUuX0KFDB4k6lz86Tkk5l18pKSkhNzcXnTt3xvjx4+Hs7IzGjRtX+XxS0qlBsrKyoKurCzs7O0RFRWHTpk04fPgwDh48KO7QflnTpk3LlOXn5wMA6tatyyurW7cur928JirvOAFAQUEBhoaG2LdvH86fP486depg2rRpIo5OsIqKiuDh4QE5OTkMHjwYgGSdy6++PU5HR0eJPJd169bF3bt3cfz4cQQFBfG+c6pyPmv11AY1TdOmTXHo0CHe+06dOsHZ2Rnh4eFwdnYWY2TC8fUH+tsf4oKCAgBfvqQliZeXF9/7hQsXomfPnkhPT6/00PHVQUZGBmbOnAkA2LNnD6Skvvx9K2nn8r/HqaCgIHHnEgCkpaUhLS0NXV1djBw5Evfv3wdQtfNJVzo1yOPHj+Hv789X9vnzZ8jJyYkpIuGqX78+VFVVkZKSwitLSUmBpqZmpYdTr84YY/D19eU7zq9NFHXq1BFXWFX25MkTjBgxApqamti/fz9UVFQk8lyWd5ySdi5jY2MxevRovrLi4mIoKytX+XxS0qlBFBUVERgYiJMnT4LL5eL+/fs4ePAgRo4cKe7QhMbGxgZbtmzB+/fv8fbtW2zfvh3Dhg0Td1gCxeFw8PDhQ6xbtw6fPn3Chw8fsGrVKpiamkJVVVXc4VVKdnY2Jk6cCGtra3h7e/P9QSRJ5/J7xylJ5xIAOnbsiBcvXmD//v0oLS1FfHw8QkJCMHLkyCqfT5pPp4a5cuUKNmzYgNTUVKioqMDV1VWiZj6Ni4uDm5sb7ty5A+DLldzatWtx4cIFlJSUwNbWFgsWLKjxk+799zjfv3+PVatW4dq1aygtLYWpqSmWLl2K+vXriznSytmzZw+8vb1Rt25dcDgcXnmfPn3g4+MjMefyR8e5bNkyiTiXXyUkJODPP/9EUlIS1NXVMXPmTJibm1f5d5OSDiGEEJGh5jVCCCEiQ0mHEEKIyFDSIYQQIjKUdAghhIgMJR1CCCEiQ0mHEEKIyFDSIRJtwIABaN++fbmv2i4uLo73Wfzxxx9IT09H+/bt8eTJkzLrBgcHw9jY+Kd1+vr68uqMjo4WRtikhquZ408QUgkeHh6ws7MTdxjVVkREBFRUVJCTk/PLdU2ZMgVOTk4wMTH59cCIRKKkQySeoqJijRyCRFQaNmwIRUVFgSQdRUVFKCoq/npQRGJR8xqp1TZv3oxJkybB1dUVhoaGvKm/d+3ahf79+0NfXx9jx47F3bt3eduUlJRg/fr16NGjB4yNjREYGAg7OzsEBwcDAJydnbF27Vq+/fy3uelH9Xt6emL58uXw9PSEvr4+evbsyTd3C5fLxZYtW9CvXz/o6+vDxcUFT58+xbt379CpUyfExsby1mWMYcCAAThx4oTAPrPg4OBymyslcaRzIniUdEitd/nyZRgZGeHYsWMwNTXF0aNHcfDgQSxfvhwhISEwNTXFuHHjkJ6eDuBLogoJCcG6deuwd+9eREVFISkpqcL7+1n9AHDixAmoq6vj1KlTmDBhAvz8/PDo0SMAgL+/Pw4cOIBFixYhJCQEqqqqmDJlCho1aoQePXrg3LlzvHru3LmDzMxMWFpaCujTAgYPHowrV67wXhs3boSMjAwmTZoksH0QyUVJh0g8b29v6Ovr873i4+N5y+vWrYspU6ZAS0sLDRs2REBAADw8PGBqagpNTU24ubnB0NAQhw4dAmMMx44dw/Tp09G3b1907NgRf/31F0pLSyscz4/q/6ply5aYNWsWWrZsyZsSOCEhAYwxHDlyBNOnT4e5uTk0NTWxbNkyDBo0CB8+fICNjQ3Cw8NRXFwMADh79iz69+8PJSWlCsc3atSoMp/XsmXLeMvl5eWhqqoKVVVVlJSUYOXKlXB3d4epqWmF90FqL7qnQyTe1KlTYWNjw1f27QyeGhoavJGC8/Ly8OrVK3h5eWHp0qW8dYqKiiAnJ4esrCxkZWWhU6dOfNurqalVKJaf1f9Vy5Yt+bZTVFRESUkJsrOzkZWVhS5duvCWKSsrw9PTEwAwaNAgLF++HNeuXYOJiQn+/vtv/PHHHxWK7Ss/Pz9oamrylYWHh2Pnzp18ZZ8/f4a7uzv09PRq/MyYRHQo6RCJp6KiglatWn13+beTa329YvH29uZLLMCXv/Dl5eUBfLmv8q0fTaRXUlJS4fp/VB9jDLKysrz/l0dRUREDBw7E33//DRkZGd7Q+pXRrFmzMp9Xo0aNyqy3ZMkS5OfnY926dXzD+xPyI9S8Rsg3vs6I+PbtW7Rq1Yr32rdvHy5fvgxFRUVoaGjg3r17vG2ys7Px8uVL3ns5OTnk5eXx3qelpVW4/p+pV68eGjVqxLu/A3yZMrh37968KYRtbGwQHR2NyMhIWFhY8BKVIO3duxcRERHYsmVLpZruCKErHUL+Y9KkSdi6dSuaNGmCzp074/Tp0wgKCsKBAwcAfHkWZcOGDWjRogVat26N9evX86YkBoDOnTvj+PHjsLW1haKiYpkZNH9W/8+MHz8e27ZtQ/PmzdGqVSts3boVSkpK6NChAwDAxMQE0tLSOHHiBHbv3i3AT+aL2NhYrF+/HqtXr0aDBg2QkZHBW0Zd08nPUNIh5D9cXFxQWFgIHx8fZGZmonXr1vDz84OBgQEAwN7eHjk5OVi+fDk+f/4MV1dX3gygADBx4kSkpKRg4sSJUFFRwaxZs/iuhH5W/89MnDgReXl58PLyQl5eHgwMDBAQEMBLbNLS0rCyskJUVBQMDQ0F+Ml8cerUKRQXF2P+/Plllj1+/Fjg+yOShWYOJUQAjI2NsWDBgmoz8sGsWbPQunVrzJ0797vrxMXFwcXFBbdv3xb4A53t27fH9u3b0b9/f4HWS2o+utIhRILcuHEDjx49QkxMTLlXIuXJysoCY0wg92by8vKQn5//y/UQyUUdCQiRIGFhYfDz88P8+fPRokWLCm1jZmaGDRs2CGT/gYGBNO4a+SFqXiOEECIydKVDCCFEZCjpEEIIERlKOoQQQkSGkg4hhBCRoaRDCCFEZCjpEEIIEZn/A31TGgqsdDHXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the aperiodic component on a linear-log scale\n",
    "plt.plot(freqs, psd_aperiodic[1, :], 'k', lw=2)\n",
    "plt.xlim(1, 30)\n",
    "plt.yscale('log')\n",
    "sns.despine()\n",
    "plt.title('Aperiodic component, chan = ' + chan[1])\n",
    "plt.xlabel('Frequency [Hz]')\n",
    "plt.ylabel('PSD log($uV^2$/Hz)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEgCAYAAACXa1X+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABLNElEQVR4nO3dd1hT59sH8G8IewkisgVFmY4CCiooiqitA/eoW7Q4axXrrKOtC8VVxb13qxUcde+9F4rKVAuKGGUIyEye9w/fnJ8poCEEMrw/18WlOefknPuQkG+e5zmDxxhjIIQQQuREQ9EFEEIIUS8ULIQQQuSKgoUQQohcUbAQQgiRKwoWQgghckXBQgghRK4oWFRIXl4eIiIi8O2336JBgwZo0aIFpk6diufPn8tl/StXrkT37t0BAJGRkfDx8QEA3LhxA87OzsjNzZVqPU+fPsX169flUhMpn/T0dBw+fFgu65o6dSrGjRsnl3VVphcvXmDatGlo0aIFGjVqhM6dO2PXrl349EwKZ2fnMn8iIyO/uI2IiAg4OzvjyJEjJeaJRCKsX78egYGB8PDwQJ8+fXDnzh257qOqoWBREXl5eRg0aBCOHj2K0NBQHD9+HCtXrkRhYSF69uyJe/fuVXgbwcHB2LRpU4XXM3r0aCQkJFR4PaT8wsPDcerUKUWXUWWio6PRo0cPCIVCrFq1CocPH8aQIUOwbNkyzJ8/n1vu8uXLJX46deoEAwMDNGnS5IvbOXz4MBwcHLB///4S87Zv345NmzZh6tSpiIqKgpeXF4YPH46XL1/KdV9ViaaiCyDSWblyJTIyMhAVFQUjIyMAgI2NDb755htMnz4dkydPxtGjR6GlpSXzNgwMDGBgYCCvkokCfE3nOzPGMHXqVLRp0wYLFy7kpteqVQu6urqYOHEi+vbtC0dHR5ibm0s89/Dhw/jnn3+waNEi2NnZfXY70dHRePHiBRYvXoxJkybh1atXsLa25ubv378fw4YNQ2BgIABg8uTJOHPmDE6dOoUhQ4bIb4dVCLVYVIBIJMK+ffswePBgLlQ+9dNPPyE5ORlXrlwBANy+fRs9e/ZEw4YN4efnh0WLFkEoFHLrWrVqFVq1agUPDw8MGjQIiYmJACS7wj4nOjoagwcPhoeHBxo0aIBevXrh/v37AICBAwfi5cuXmDNnDgYOHAgAePv2LaZMmYJmzZrB09MTP/30E968eQMASElJgbOzM1avXg1vb2/88MMP6NChA5YvXy6xzSlTpmDy5Mml1vPq1SuMHj0anp6eaN68OebPn4/i4mIAQE5ODubOnYuWLVuiUaNGGDZsGJKSkrjnBgQEYN++fejXrx8aNmyIHj164Pnz51i4cCEaN24MPz8//PXXXxLLb9myBX369EHDhg3Ru3dvPHz4kJtfWFiIiIgIBAQEoEGDBvj+++/x4MEDbv7AgQOxYsUKjBo1Cg0bNoS/vz/27dvHzS8qKsLixYvh6+sLLy+vEvV+7vkrV65EVFQUTpw4AWdn5y++jgBQUFCA+fPnw9fXF56enhg9ejTS0tK4+Xl5eZg+fTo8PT3RtGlTideluLgYy5YtQ0BAANzd3dG8eXPMmzePe6+tXLkSY8eORVhYGLy9vdG4cWPMmTOHm/9fAwcOLLO7qjT3799HYmIiQkJCSsz79ttvsXXrVtja2paYl5ycjNmzZyMoKAhdunT54u/o0KFDqFevHtq3bw8DA4MSXWczZ85EUFBQiedlZ2d/cd1qixGll5CQwJycnNiDBw/KXKZ9+/Zs6dKlrLi4mHl7e7Pw8HCWnJzMrl69yho3bsz27t3LGGPsjz/+YD4+PuzEiRPs2bNnLDQ0lAUEBLDi4mK2YsUK1q1bN8YYY/v372fe3t6MMcauX7/OnJycWE5ODsvOzmbe3t5swYIF7MWLF+zx48ds0KBBrHPnzowxxjIyMljLli3Z2rVrWUZGBisqKmIdO3Zk/fv3Zw8fPmTR0dGsT58+rGfPnkwkErHk5GTm5OTE+vTpw549e8bi4uLYmjVrWLt27bh9y8/PZ56enuzSpUsl9rugoIC1a9eOBQcHs8ePH7Nbt24xf39/tnz5csYYY8HBwaxTp07s1q1b7OnTp2zUqFGsVatW7MOHD4wxxlq3bs18fHzYmTNnWFxcHGvXrh1r0qQJmzt3LktMTGQLFy5k7u7u7N27d9zyjRo1Yjt37mQJCQlsypQprEmTJiwjI4MxxtjMmTNZy5Yt2fnz51lCQgKbOXMm8/DwYGlpaYwxxgYMGMDq16/Pdu3axf799182Z84c5u7uzgQCAWOMscWLF7POnTuzmzdvsoSEBDZnzhzWvHlz9v79+y8+Pycnh/30009sxIgR7M2bN1K9tyZPnswCAgLYlStXWEJCAhs8eDDr168fY4yxKVOmMCcnJ7Z06VL277//sr179zInJyd2+fJlxhhja9asYf7+/uzGjRssOTmZRUVFMVdXV3by5EnGGGMrVqxg7u7ubMaMGSwpKYlFRUUxFxcXbv5/ZWRksDdv3pT6U5pdu3YxDw8PqfZTrKioiPXu3Zu1adOGZWdnS7V8s2bNuPfTxIkTWUBAABOJRGU+59y5c8zJyYndunWrXLWpEwoWFXDnzh3m5OTEXrx4UeYyffr0YTNnzmQZGRnM2dmZbdmyhXvzP3jwgKWkpDCRSMSaNm3Ktm3bxj0vKyuLLViwgL17906qYBEIBGz9+vWsqKiIW8fx48eZi4sL97h169Zsx44djDHGzp49y9zd3dnr16+5+ampqczNzY1dvnyZC5YjR45w81NSUpizszN79OgRY4yxY8eOMV9fX1ZcXFxiv8+dOyfxwc8YY+fPn2e7d+9msbGxzMnJiUVHR3PzcnNzmbe3N/vrr7+4WufMmcPNX7RoEfP29ua29e7dO+bk5MTu3r3LLT9t2jRu+YKCAtasWTO2a9culpWVxVxdXdmxY8e4+UKhkHXo0IEtXbqUMfYxGIYNG8bNz87OZk5OTuzSpUssLy+P1a9fn925c0diH9u1a8d27tz5xecz9jEMfvzxxxK/p9K8f/+eubm5sdOnT3PTXrx4wcLDw1lBQQGbMmUK69Spk8RzAgMD2fr16xljjJ06dYpdvXpVYn6nTp1YREQEY+xjsHh5ebGCggJufteuXdmSJUukqu9L1qxZw1q0aFGu5yxbtoy5ubl99kvap86fP8+cnJzY48ePGWMf99nJyYldu3at1OXj4+OZt7c3Cw0NLVdd6obGWFRAtWrVAHzs1ilLdnY2TE1NYWJigkGDBmHBggXYsGEDWrZsiQ4dOqBhw4ZIT09Heno6GjRowD3P2NgYU6dOlbqWGjVqoHfv3ti9ezeePn2K58+f4/HjxxCJRKUun5CQAGtra1hYWHDTLC0tYWNjg/j4eNjb2wOARJeFjY0NGjdujKNHj8Ld3R1HjhxBhw4dwOfzS12/jY0Nqlevzk3z9/cHAG7MqX79+tw8fX19uLm5IT4+npv26bZ1dXVhbW3NbUtHRwfAxy4uscaNG3P/19bWhouLC+Lj4/H8+XMIhUJ4eHhw8zU0NODh4SGxPQcHB+7/hoaGAD52K/37778oLCxEcHAweDwet0xBQYFEd1hZzy+vZ8+eobi4WOL9UKtWLfz8888Sjz9lZGSEgoICAEBgYCBu3ryJ8PBwPHv2DHFxcUhOTka7du245a2srKCtrS1Rb1FRUan1DB8+vMyjqUo7OMXU1LRc3U23b9/GunXrMH78eDRs2FBi3qxZsySOptuwYQMaN26MQ4cOwc7ODq6urgCAFi1aQF9fH/v370fTpk0l1vHo0SP88MMPqFevHubNmyd1XeqIgkUFODg4oHr16rhz5w7c3NxKzE9PT8ezZ8/QqFEjAMD06dPRr18/nDt3DufPn0dISAjGjBmDwYMHA6jYAG9aWhp69eqF2rVrc6GVkZEh8WH0KfEH83+JRCKJMNLV1ZWYHxQUhLVr12LUqFG4cOECdu/eXep6PnewgrTb1tSU/DP49EO9NP8NOKFQCD6fL/EB+rntlVYzY4wbe9i8eTPMzMwk5osD5HPPLy9xvZ97bmlhLl4+IiIC27ZtQ48ePdCuXTtMnDgREydOlFi2PAeTzJs3D/n5+VIv36BBA3z48AHx8fGoV6+exDyhUIgRI0Zg6NCh8PX1xfv37zFp0iQ0adIEP/zwQ4l1/fTTTxg2bBj32MLCArm5uThz5gzy8/Ml/u6EQiFOnjyJ2bNnc6/LrVu3MGLECDRq1AirV68u8X7+2tDgvQrg8/no27cvNm/ejMzMzBLzly9fDisrK7Ro0QJv377Fb7/9hho1amDo0KHYtm0bQkJCcOTIERgZGcHMzAyPHz/mnpufnw9fX19ER0dLVcuRI0egqamJrVu3YtiwYfDz88Pr168BlP4B5ejoiFevXkkMCL9+/RqpqalwdHQsczvffvstBAIBtmzZAhsbG4lWx6ccHBzw8uVLZGRkcNOioqLQq1cvODo6oqioSGJw/cOHD3j69Cnq1Kkj1f6WJiYmhvt/fn4+YmNj4ezsDHt7e2hpaeHu3bvcfMYY7t+/L9X2atWqBU1NTaSnp8Pe3h729vaws7PDH3/8IXEAwOd8KRQ/ZWtrCz6fL/F+SE5ORtOmTSVer7Js2rQJU6ZMwdSpU9G1a1fY2tri1atXMn9xsbCw4Pb7vz+lcXNzg5OTEzZu3Fhi3pEjR3Dp0iWupTxr1ix8+PAB4eHh0NAo+bFnZmYmsT1dXV2cOnUKeXl5WLduHQ4cOMD9LFmyBPn5+fjnn38AfDxvKyQkBE2aNMG6deugp6cn0/6rE2qxqIhRo0bh1q1b+P777xEaGgp3d3e8ffsWO3fuxJkzZ7BhwwZoaWmhWrVqOH36NPLz8xESEoK8vDxcvnyZa/oPGTIEa9asga2tLezt7bF69WoYGhrCxcUFFy5c+GIdFhYWePv2Lc6fP4969erh5s2bWLNmDYCP3UU6OjowMDBAYmIi3r17h+bNm8PZ2RkTJ07E1KlTwRjDggUL4ODggGbNmnFHh/2XsbExWrdujQ0bNmDkyJFl1uPn5wd7e3tMmzYNoaGhyMrKwsqVK9GtWzc4ODigXbt2+OWXXzB79mwYGRkhIiICmpqa6Nixowyvwkd79+5FgwYN4ObmhtWrV0NLSwvfffcd9PT0MGDAAISFhUFfXx+2trbYsWMHXr58id69e39xvQYGBvj+++8xb948aGlpwcHBAZs2bcK5c+cwYcIEqWrT19dHXFwcUlJSYGtri9zcXHz48KHE4bbAx1ZQz549ERYWBkNDQ5iYmCAsLAz16tWT6Losi4WFBS5cuIAmTZogJycHK1euRFZWlkS3YWX77bffuK7D/v37w9DQEJcuXcKyZcsQHByMunXrIioqCseOHcPcuXOhoaEBgUAgsQ59ff1SD7M/dOgQPDw8uK5VMScnJ6xbtw779+9H3759MW3aNNSsWROzZs1CVlYWt5yenp5ES/NrQsGiIrS1tbF582Zs3boVy5cvR0pKCoyNjeHr64v9+/dz/e5aWlpYv3495s+fj+7du0NTUxMBAQGYMWMGgI8nQebm5mLGjBnIzc2Fp6cn1q1bV2Y3zn999913uHfvHqZOnYqioiLUq1cPc+fORWhoKB49egQvLy8MGjQIYWFhuHfvHg4cOIDVq1dj3rx5GDhwIPh8Plq2bIk//vjji9vs1KkTTpw4gc6dO5e5DJ/Px5o1azBnzhz06tULRkZG6NatG0aPHg0AmD9/PhYsWIBRo0ahuLgYTZo0wc6dO2FiYiLV/pamZ8+e2LJlCxITE/HNN99g69at3AdIaGgoeDwepk2bhtzcXDRo0ADbt2+XGBf5nMmTJ0NTUxPTp09HTk4OXF1dsXHjxi+eayHWrVs3nD59Gh07dsTp06fx559/IiIiArGxsaUuP3XqVISFhWHs2LEQCoXw9fXFzJkzpdpWWFgYfvvtN3Tu3BnVq1dHYGAgevbsKdGiq2yenp7YtWsX1qxZg5EjRyI3Nxf29vaYOnUqevXqBeDjlSMAcH8D/zV27Fj8+OOPEtPevHmD69evlzlWMmDAAMyaNQtPnz7lWnwBAQESywQHB2PKlCkV2j9VxWMV6XAnpBLt2rUL//zzD/bs2aPoUjgBAQEIDg7GgAEDFF2K1Lp37y7VZUsIkRdqsRClExcXh9jYWKxdu7bMgwKIdI4ePcod1EFIVaFgIUrn6dOnmDlzJtq1a/fZbjDyZYGBgejQoYOiyyBfGeoKI4QQIld0uDEhhBC5UvtgKS4uRkpKikxnJhNCCCk/tQ+W169fo02bNtxJfIQQQiqX2gcLIYSQqkXBQgghRK4oWAghhMgVBQshhBC5omAhhBAiVxQshBBC5IqChRBCiFxRsBCl9OzZs3LdTZAQojwoWIjSefr0Kfz8/DB9+nRFl0IIkQEFC1E68fHxAIDExEQFV0IIkQUFC1E6mZmZAID3798rthBCiEwoWIjSoWAhRLVRsBClIw6W7OxsxRZCCJEJBQtROuJgyc3NpdsdEKKCKFiI0hEHC0CtFkJUEQULUToULISoNgoWonQ+DRYawCdE9VCwEKVDwUKIaqNgIUqHusIIUW0ULESp5OXlSVwjjFoshKgeChaiVD5trQAULISoIgoWolSysrIkHlOwEKJ6KFiIUvlvi4XGWAhRPRQsRKlQsBCi+jQVXYBYdHQ0wsLCEBcXB0NDQ/To0QNjx45FQUEBvLy8oK2tzS3r4eGBzZs3K7BaUlnEwaKrq4v8/PwSXWOEEOWnFMGSm5uLESNGYMyYMdixYweSk5MxbNgwVK9eHe7u7jAzM8PFixcVXSapAuJgqVWrFuLi4qjFQogKUopgSU1NhaenJwYMGAAAcHBwQNu2bXH37l0AgKurqyLLI1UoIyMDAGBnZ0fBQoiKUooxlrp162LVqlXc48LCQly8eBEuLi6IiYmBQCBA586d0bx5c4wbNw5paWkKrJZUJnGLxd7eHkDJo8QIIcpPKYLlU4WFhQgNDYW2tjb69+8PfX19eHl5Ydu2bTh27Bh0dHQwatQoRZdJKok4WOzs7ADQ4D0hqkgpusLEBAIBfvzxRwDAli1boK+vjxkzZkgsM23aNDRr1gwpKSmwtbVVRJmkEn06xgJQsBCiipSmxRIXF4cePXrAwcEB27dvh6mpKRhjWLZsGZKSkrjlCgsLAQA6OjqKKpVUInGwWFpaQktLC/n5+SgoKFBsUYSQclGKYMnIyEBwcDA6deqEsLAw7tBiHo+HmJgYLFq0CNnZ2cjKysLcuXPh7+8Pc3NzBVdNKoM4WExMTGBkZASAWi2EqBqlCJYDBw5AIBBg9+7d8PDw4H7GjRuHhQsXQk9PD4GBgWjTpg10dHQQHh6u6JJJJfk0WIyNjQHQZV0IUTVKMcYydOhQDB06tMz5y5Ytq8JqiKIUFRUhJycHPB4PxsbG1GIhREUpRYuFEOB/LZNq1apBQ0ODa7HQIceEqBYKFqI0xCdHmpiYAAAXLNRiIUS1ULAQpSEeXzE1NQUA6gojREVRsBCl8enAPfC/YKHBe0JUCwULURr/DZZq1aoBoGAhRNVQsBClQS0WQtQDBQtRGv8NFhq8J0Q1UbAQpVFWsFCLhRDVQsFClAZ1hRGiHihYiNIQB4t40J66wghRTRQsRGmU1WKhYCFEtVCwEKUhPvNefIKkuOVCl3QhRLVQsBClIQ6Q0losjDFFlUUIKScKFqIURCIRFyziloq2tjZ0dXVRXFyM/Px8RZZHCCkHChaiFLKzsyESiWBoaAgtLS1uurjVQt1hhKgOChaiFP47cC9GR4YRonooWIhSEAeHuIUiRidJEqJ6KFiIUsjNzQUAGBgYSEynQ44JUT0ULEQpiIPF0NBQYjrdRZIQ1UPBQpSCOFj09fUlptMYCyGqR1PWJ4pEIiQlJeHdu3fg8/moUaMG7O3twePx5Fkf+UpQVxgh6qPcwXL79m1s374dV65c4T4MAIDH48HY2BgtWrRAv3794OnpKddCiXr7UrBQVxghqkPqYHnx4gVmzZqFV69eITAwEBEREXB0dISJiQlEIhEyMjLw9OlT3Lp1C6GhobCzs8Pvv/+O2rVrV2b9RE2UFSzikyWpxUKI6pA6WCZNmoQxY8bA39+/1PlWVlawsrJC69atMWnSJJw5cwaTJk3C33//LdX6o6OjERYWhri4OBgaGqJHjx4YO3YsioqKMG/ePBw/fhwA0LVrV0yePBl8Pl/a0okKKCtYzM3NAQCvXr2q8poIIbKROlj++usvqcdPeDweAgMD0aZNG6mWz83NxYgRIzBmzBjs2LEDycnJGDZsGKpXr47Xr18jISEBJ06cQH5+PkaOHIlNmzYhJCRE2tKJCvjw4QOAksFSt25dAEB8fHyV10QIkY3UR4V9GirTp09HbGxsiWXS09MlwkTaIEpNTYWnpycGDBgAPp8PBwcHtG3bFnfv3kVUVBRCQkJgYmICS0tLjBw5Evv375e2bKIiymqx1KlTBxoaGnjx4gUKCgoUURohpJxkOtw4MjISffv2xT///CMxXSQSydRlUbduXaxatYp7XFhYiIsXL8LFxQUCgQCOjo7cvDp16uDFixcoLCyUpXSipHJycgCUPNxYR0cHtWrVgkgkwrNnzxRRGiGknGQ+j2X8+PGYOXMmFixYAJFIJLeCCgsLERoaCm1tbXTo0AEAoKenx83X09MDY4yudqtmyuoKA4B69eoBABISEqq0JkKIbGQOlk6dOmH37t04deoUhgwZgvT09AqfwyIQCDBo0CC8ffsWW7Zs4c7C/jRE8vLyAJT8ZktUW1ldYQCNsxCiamQKFnGAuLq6ckd9de/eHXFxcTIXEhcXhx49esDBwQHbt2+HqakpqlWrBnNzcyQlJXHLJSUlwcHBAZqaMp/bSZRQWZd0AajFQoiqkSlYPr2bX/Xq1bFlyxYEBgZixIgRMhWRkZGB4OBgdOrUCWFhYdDW1ubmBQUFYdWqVXj37h3S0tKwdu1adO3aVabtEOVV1iVdAGqxEKJqZPraP3bsWIkPAD6fjxkzZsDd3R2RkZHlXt+BAwcgEAiwe/du7Nmzh5veokULhIeHY+HChQgKCkJxcTG6dOlChxqroc91hYlbLImJiRCJRNDQoEvcEaLMeEzKm4mXpxtC/A1TGaSkpKBNmzY4c+YMbG1tFV0OKUO9evXw4cMHPH36tMQ9WQDA09MTaWlpuHbtGmrVqqWACgkh0pK6xdKpUyfweDwwxiT+BSS7xng8Hp48eSL/SonaEolE3FFhZR2U4ejoiLS0NCQkJFCwEKLkpA6WM2fOSDxmjKFz585Yv349rK2t5V4Y+XqIj/TT09Mr81I99erVw9WrVxEfH4+AgICqLI8QUk5SB4uNjU2p0y0tLcucR4g0Pje+IkZHhhGiOmgUlCicNMFCR4YRojooWIjCfe5QY7FPg0XK400IIQpCwUIUTpoWi6WlJQwNDZGZmYn09PSqKo0QIgOpx1h27dpVYppQKMTBgwdhamoqMb1///4Vr4x8NT531r0Yj8dDvXr1cO/ePcTHx8PMzKyqyiOElJPUwbJp06YS02rUqFHihEgej0fBQspFmq4w4GN3mDhYmjZtWhWlEUJkIHWwnD17tjLrIF+xz13Z+FPi81dSU1MrvSZCiOykDpa2bdvCz88Pvr6+aNq06We7LQgpD/G9WL4ULOIWjTiICCHKSerB+/Xr16Nu3brYv38/AgIC0LdvX0RERODu3btyvR8L+fpIM3gPULAQoiqkbrHUrl0btWvXRv/+/SEUCnHv3j1cvXoVCxcuxPPnz9G4cWP4+vrC19cX9vb2lVkzUTMULISoF5mubszn89G4cWM0btwY48aNQ3Z2Nq5du4YrV65g69atOHnypLzrJGpM2jEW8XwKFkKUm1zulmVkZIR27dqhXbt28lgd+cpQi4UQ9SLVGEtiYiKuXbvGfQCInTt3rlKKIl8XaQ83Fs//7/uQEKJcvhgs27dvx5gxY7Bjxw506tRJoptrxYoVlVoc+TpI22LR09MD8L+rIRNClNMXu8L27duHyMhI6OvrIyUlBePGjcPLly8xdOhQumYTkQtpzrwHaIyFEFXxxWBhjHFdELa2ttixYwfGjRuH169fU7AQuaAxFkLUyxe7wqpXry5xR0gDAwOsW7cOGRkZiIuLq9TiyNehvGMsFCyEKLcvtlgWLVpU4q5+mpqaWLRoEfr06VNphZGvh7SHG38aLJ/eGpsQoly+2GKxtLSEubk5AODgwYMS87y8vCqnKvJVkfaSLnw+H7q6umCMIT8/vypKI4TIoFz3Y5k+fToGDRqExMTEyqqHfGVEIhHXYvlSVxjwvyPDqDuMEOVVrmDZv38/hEIhunTpgvDwcDrsk1SY+D2kp6dXosu1NDTOQojyK1ewuLi4YNeuXZg3bx4OHTqE7777DidOnJBrQdHR0fD29uYe5+fnw93dHR4eHtxPcHCwXLdJFEfaI8LE6CRJQpSfTJd06dKlCwIDAxEREYGJEydi3759mDVrFne/DFkdPXoUs2bNglAo5KY9ffoUZmZmuHjxYoXWTZRTeYOFzmUhRPnJfM97AwMDTJkyBQcPHoRIJEKnTp2wYsUKFBYWyrS+ZcuWYePGjRg9erTE9JiYGLi6uspaJlFy0h5qLEZjLIQoP5mDRczR0RGbN2/GkiVLcODAAXTs2FGm9fTr1w+RkZFwd3eXmB4TEwOBQIDOnTujefPmGDduHNLS0ipaNlESsnaFUbAQorxkvrrxy5cvERcXh7i4OMTGxiIuLg5paWky3/TLwsKi1On6+vrw8vLCqFGjwOfzMXfuXIwaNQqRkZGylk6UiLSXcxGjYCFE+ZUrWGbNmoW4uDjEx8dzf9h2dnZwcXHBd999BxcXF7i4uMi1wBkzZkg8njZtGpo1a4aUlBTY2trKdVuk6pW3K4zGWAhRfuUKltjYWDg7OyMoKAguLi5wdnaWugtDFowxLF++HF26dEGdOnUAgBvD0dHRqbTtkqoj7Vn3YtRiIUT5lStY/vrrr8qqo1Q8Hg8xMTGIjY1FeHg4RCIR5s6dC39/f+5qAES10RgLIepHpsH7AwcOID4+XubxlPJYuHAh9PT0EBgYiDZt2kBHRwfh4eGVvl1SNaS9nIuY+KgwOo+FEOUl0+B9WFgYsrKyoK2tjXr16sHNzQ2urq5wc3ODs7MzdHV1ZS7Ix8cH9+7d4x6bmZlh2bJlMq+PKDdZz2Ohqz4QorxkCpbr168jNTUVMTExePLkCR48eIDIyEgIhUJoaGggJiZG3nUSNUVjLISoH5kPN7aysoKVlRUCAwMBAK9fv8aECRPQpUsXuRVH1B+NsRCifip8gqSYpaUlFixYgAsXLshrleQrUN7DjelaYYQoP5mCJT09vdTpDg4OePz4cYUKIl8XarEQon5k6gpr3rw5LCwsuAF7Nzc31KlTB9HR0SgqKpJ3jUSNyXrmPQ3eE6K8ZAqWPXv24MmTJ3jy5AnOnz+PjRs3orCwEHw+H5MnT5Z3jUSNUYuFEPUjU7CI74siJhQK8fLlS1SrVg3VqlWTW3FE/ck6xkLBQojyknqM5f3792XO4/P5qFWrVolQ+dxzCAHKf7ixeDkavCdEeUkdLAMGDMD69euRnZ39xWXT09OxevVq9O/fv0LFEfVHXWGEqB+pu8L27NmD5cuXo2XLlvD09ESLFi1Qr149mJqaQiQSISMjA0+fPsWtW7dw69YtdO/eHXv27KnM2omKE4lE5e4KE1/VIT8/HyKRCBoacjtinhAiJ1IHi4GBAX755ReEhIRg7969OHLkCB4/fszdRlhTUxNubm7w9/fHnDlzyry/CiFi4q5SfX198Pl8qZ6joaEBPT095OXlIS8vr1Kvrk0IkU25B+/Nzc0xZswYjBkzBiKRCJmZmeDxeDA1Na2M+ogai4uLAwDulgjSMjAwQF5eHnJzcylYCFFCMl/SBfj47bF69eryqoV8ZcQn07q5uZXreTTOQohykylYXFxcwOPxSp2npaUFCwsLdOjQAePGjZO6i4N8fShYCFFPMgXL77//jj/++AM//vgjvvnmGwDAw4cPsWLFCvTu3Rv16tVDREQEAGDChAlyK5aoF1mDRXxPFgoWQpSTTMGyadMmzJs3D61ateKmubi4wNLSEvPmzcPx48dhbm6O0NBQChZSKqFQiCdPngAof7DQfe8JUW4yHav55s0b2NnZlZhuaWmJ1NRU7v9ZWVkVq46orWfPniE/Px9WVlblPvCDusIIUW4yBYunpyfCw8ORmZnJTcvMzMSSJUvg6ekJADh16hQcHBzkUSNRQ7J2gwEULIQoO5m6wubMmYOQkBC0bNkSNjY2EIlESE1NRZ06dRAREYGLFy9i6dKlWLFihbzrJWqCgoUQ9SVTsFhbW+PgwYO4fv06YmNjoampiXr16qFZs2YAPv7hnz9/ng5FJmWSR7DQ9cIIUU4yn8fC5/OhoaEBTU1NiEQiFBUVobi4GJqamhQo5IuoxUKI+pIpWAQCAcaMGYPHjx9LdIXVrl0bW7duhZmZmbzrJGokIyMDqamp0NXVRe3atcv9fLrZFyHKTabB+3nz5oHP5+PMmTM4ceIETp06hdOnT8PY2BgLFiyQd41EzYhbK66urjKdQEstFkKUm0zBcunSJcyYMUPiQpOWlpaYOnUqLl68WKGCoqOj4e3tzT0uLCzE7Nmz4ePjAx8fHyxYsIC78CVRTRXpBgPoniyEKDuZgkVHR6fUy5VraGhU6EP/6NGjCA4ORlFRETdt5cqVSEhIwIkTJ3Dw4EHcuHEDmzZtknkbRPEqGix05j0hyk2mYGnRogXmz5+Pt2/fctMEAgHCwsLg5+cnUyHLli3Dxo0bMXr0aInpUVFRCAkJgYmJCSwtLTFy5Ejs379fpm0Q5SAOFnd3d5meT11hhCg3mYJl8uTJyMjIQOvWrdG+fXu0b98eAQEByMvLwy+//CJTIf369UNkZKTEh8379+8hEAjg6OjITatTpw5evHiBwsJCmbZDFCs9PR1PnjwBn8+Hq6urTOugwXtClJtMR4WZmZnh4MGDuHjxIhITE6GjowNHR0c0b95c5kJKuzGY+BupuOtD/H/GGPLz86GtrS3z9ohiHDt2DEKhEK1atYKhoaFM66AWCyHKTepgSUhIKDHNzs5O4pph4mXq1q0rh9L+Fyj5+fncNPG3VGlvZUuUy6FDhwAAQUFBMq+DBu8JUW5SB0unTp3A4/HAGOOm/feeLIwx8Hg87qq1FVWtWjWYm5sjKSkJNjY2AICkpCQ4ODhAU7NC9ygjCiAQCHD16lVoaWmhffv2Mq+HWiyEKDepP53PnDlTmXWUKSgoCKtWrYKbmxuKi4uxdu1adO3aVSG1kIo5cuQIRCIRWrduDRMTE5nXQ8FCiHKTOljELYaq9tNPP2HhwoUICgpCcXExunTpgpCQEIXUQirm8OHDACrWDQb8r4uUBu8JUU489mnflhpKSUlBmzZtcObMGdja2iq6nK/W69ev0bhxY2hra+PBgwcwMjKSeV2MMdSqVQsikQgvXrygblFClIxMhxsTUl5HjhwBYwytW7euUKgAH8f2qDuMEOVFwUKqxLFjxwBUvBtMjIKFEOVV7j6E9PR0nD9/HvHx8cjNzYWhoSFcXV3h7+8PY2PjyqiRqIFnz54BALy8vOSyPgoWQpRXuYJl48aNiIiIAI/Hg42NDYyNjZGTk4MdO3ZAU1MTEyZMwKBBgyqrVqKiRCIRd/kfc3NzuayTgoUQ5SV1sOzbtw8RERGYNGkSunfvLnE2fH5+PiIjI7F48WJYWFhU6BwFon4yMjJQXFwMExMT6OjoyGWdFCyEKC+pg2Xnzp34+eef0b9//xLzdHV10a9fPwiFQmzfvp2ChUh48+YNAPm1VgAKFkKUmdSD9y9evPjilYv9/PyQmJhY4aKIehEIBAAoWAj5WkgdLPn5+V+8aKCRkRGysrIqXBRRL+IWS2kXGpWVOFhycnLktk5CiHyU63Dj/14bjBBpVEaLpVatWgAgt+vSEULkp1xHhUVGRn72qsJ0tVlSGnGLpWbNmnJbp4+PDwDg+vXrclsnIUQ+pA4Wa2tr7Nmz54vLWVlZVaggon4qo8XSuHFjaGpq4tGjR3j//j2dQ0WIEpE6WM6ePVuZdRA1lpaWBkC+LRZ9fX00atQId+7cwa1bt9CmTRu5rZsQUjF0SRdS6cQtFnkGCwA0a9YMAHWHEaJsyhUsOTk5OHDgAHckDmMMGzZswKBBgzB+/HhER0dXSpFEtVVWsIjHWa5duybX9RJCKkbqYPn333/x3Xff4ffff0dGRgYAYN68eVi6dClq1KgBY2NjDB48GHfv3q20YonqKSgoQGZmJjQ1NSt0c6/SNGnSBBoaGnj48CEdOEKIEpE6WP744w+4u7vj6tWrsLOzw7t377Bnzx58++23WLp0KX7//XeMHTsWERERlVkvUTHia4TVqFEDGhry7Xk1MjJCgwYNUFxcjDt37sh13YQQ2Un9l37t2jWMGTMGurq6AIBLly5BJBJJ3CbY19cXDx48kHuRRHVVxqHGn2ratCkA6g4jRJlIHSzZ2dkwMzPjHt+4cQOamprw9vbmphkYGEAkEsm3QqLSKuM6YZ8SBwsN4BOiPKQOFhsbGyQlJQEAhEIhLl26hCZNmkhc5fj69et0+18ioTIu5/Ipb29v8Hg83L9/H3l5eZWyDUJI+UgdLN27d8ecOXNw+PBhzJw5E2/fvsX333/Pzb99+zaWL1+Ojh07VkqhRDVVxsmRnzIxMYGrqysKCwtx//79StkGIaR8pA6WYcOGoU2bNpg7dy7OnTuH0NBQtG3bFgAwd+5cDBgwAF5eXhg2bFilFUtUT2WPsQCAh4cHAODhw4eVtg1CiPSkPvOez+dj8uTJmDx5col5vXr1Qo8ePeDq6irX4ojqq+wWCwC4u7sDAGJiYiptG4QQ6UkdLB8+fMCCBQtw4sQJaGtrIzAwED///DMMDQ3h7OxcmTUSFVYVLZb69esDoGAhRFlI3RW2cuVKnDt3DsOHD8fQoUNx8eJFzJgxozJrk7B37164u7vDw8OD+4mKiqqy7RPZVPZRYQDg6uoKHo+H+Ph45OfnV9p2CCHSkbrFcuzYMYSHh3PXZ2rSpAn69++PoqIiaGlpVVqBYo8ePcLw4cMxYcKESt8WkQ/GWKVdzuVT+vr6cHR0REJCAuLi4tCwYcNK2xYh5MukbrEIBAI4Ojpyjxs0aADGGN69e1cphf1XTEwMjeGomPfv36OgoACGhoafvY+PPFB3GCHKQ+pgEQqF4PP53GMejwctLS0UFRVVSmGfKioqQlxcHCIjI+Hn54e2bdti/fr1YIxV+raJ7Kpi4F5MPID/6NGjSt8WIeTzynUHSUVJT09Hw4YN0b17d0RERCAhIQGjR4+Gnp4eBg4cqOjySBmqYuBejFoshCiPCt2aWCgU4uDBgzA1NZVYrn///vKp7v9ZWFhg165d3GM3NzcMHDgQJ0+epGBRYoposTx+/BgikUjuF7wkhEivQrcmrlGjBiIjIyWm8Xg8uQdLbGwsTp06hbFjx3LTCgoKoK2tLdftEPmq7Mu5fMrMzAyWlpZ4/fo1nj9/jjp16lT6NgkhpVOJWxMbGBhg/fr1sLKyQrdu3fDo0SPs3LkTM2fOVFhN5Muq4lDjT7m7u+P169d49OgRBQshCqQS/QW2trZYvXo1du3aBS8vL4wfPx5jxozBd999p+jSyGdU5RgLQOMshCgLlRi8BwA/Pz/4+fkpugwipXfv3uHcuXMAADs7uyrZJl3ahRDloBItFqJ6pk+fjnfv3qF58+bcPVMqG7VYCFEOFCxE7g4dOoR//vkHBgYGWLp0aZUdoWVnZwcjIyO8efMGL1++rJJtEkJKomAhciUQCDB9+nQAwIwZM6qsGwwANDQ00KJFCwDAzp07q2y7hBBJFCxErlauXImMjAy0aNFCIecYhYSEAAC2b9+O3NzcKt8+IYSChciZ+LD0iRMngsfjVfn2mzRpAi8vL2RmZuLPP/+s8u0TQihYiBwlJyfj2bNnMDY25u7qqAijRo0CAGzYsAHFxcUKq4OQrxUFC5GbixcvAgB8fX2hqam4I9nbtWuH2rVrIzk5GUeOHFFYHYR8rShYiNxcuHABANCyZUuF1sHn87mxljVr1tBVsAmpYhQsRC6EQiGuXLkCQPHBAgC9evWCqakpHj58SOe1EFLFKFiIXERHRyMzMxP29vZwcHBQdDnQ09NDUFAQANAtrAmpYhQsRC6UpRvsU926dQMAHDhwACKRSMHVEPL1oGAhUmOMYd++fbh582aJeZcuXQKgXMHi5eUFW1tbvH79Gjdu3FB0OYR8NShYiNROnDiB8ePHo1u3bhg6dCgSEhIAADk5Obh9+zY0NDTg6+ur4Cr/R0NDA127dgVA3WGEVCUKFiK1/fv3c/8/efIkAgIC4Ovri86dO6O4uBjffPMNqlWrpsAKSxIHy5EjR1BYWKjYYgj5SlCwEKlkZWXhzJkz4PF4OH78OAYMGAAAeP78OeLi4gAA3377rSJLLJWrqytcXFyQmZmJ8+fPK7ocQr4KFCxEKseOHUNBQQGaN2+OBg0aYOHChYiJicHFixdx8OBB/P3339y5I8pG3Gr5+++/FVsIIV8JChYilcjISABA9+7duWlGRkZwdHRE48aN0axZM2hpaSmqvM/q2rUrNDQ0cOTIEUycOBF5eXmKLokQtUbBQr4oNTUVV69ehba2tkreDtrOzg5LliyBrq4u/vzzT3Tq1Alnz57Fmzdv6Kx8QioBBQv5okOHDoExhsDAQKUbnJdW7969cfjwYdSpUwdPnz7FwIED4eHhgW+++UbioARCSMVRsAAoKiqib66fIT5UVzxWoarc3Nxw/PhxjB07Ft7e3jA2Nsbbt28xefJkpKSkKLo8QtTGVxMseXl56NKlC1q3bo0FCxbg9u3b2LdvH/r164c6depg9uzZii5RKSUkJODhw4cwMjJCmzZtFF1OhRkYGGDatGmIiorC48eP0aVLF+Tn5+PXX39VdGmEqA0eU/Ov6ikpKWjTpg3atWvHDUCXRkNDA+fOnUPdunWrsDrlt2TJEixduhR9+vTB0qVLFV2O3KWmpsLf3x+5ubnYuXMnWrdureiSCFF5X02LJSoqClpaWli+fDmCg4NRu3Zt+Pj4YMGCBejRowdEIhGWLFlS6nMZY8jKyqriihWPMYYDBw4AALp06aLYYiqJlZUVJkyYAACYMWMGCgoKFFwRIWqAqYinT5+yPn36sEaNGrH27duz8+fPS/W85ORk5uTkxGxsbNgff/xR6jIpKSnMwcGBWVtbs5iYmBLz58+fz6ytrVlwcDCLj4+v0H6okujoaGZtbc0aNmzIioqKFF1OpSkoKGD+/v7M2tqaffvtt2zTpk1MIBCUWC4jI4M9e/as6gskRMWoRIulsLAQo0aNQvv27XHr1i1MmjQJ48ePx8uXL6Veh7u7O0aPHl3qPBsbGwwcOBAAEB4eLjEvLS0NGzZsAAAcP34crVu3xrRp076Kb7YHDx4EAHTq1Emhd4SsbNra2ggPD4exsTGio6Mxc+ZMeHp6omfPnli3bh1OnTqFMWPGwNPTE35+fiXuSpmeng6BQFDm+gUCAaZNm4bLly9X9q4QohRUYozl0qVLmDJlCq5cuQIejwcAGDFiBDw8PDBy5MjPPlc8xrJp0yb4+fmVuZxAIECzZs2Ql5eHw4cPw9PTEwDw22+/Yf369QgICICNjQ12794NoVCICRMm4Oeff5bfTioZkUgEHx8fvHr1CgcOHECTJk0UXVKly8vLw8mTJxEZGYnz58+juLi41OUMDAzwzz//wMnJCbdu3cKQIUPw/v17tGrVCv369UNgYCB3smheXh569OiBBw8eQFdXF1FRUWjYsGFV7pZCMMbw5MkT3LhxA9evX8fr169RvXp11KhRAzVq1ICVlRUsLS1Rp04d1K5dG3w+v1LqKCwsxP3795GdnQ1jY2MYGxvDyMgIRkZGMDAwgIZG2d+tRSIRPnz4gOzsbGRnZyMzMxNpaWkQCAR4+/YtMjMzkZmZCT6fD2tra1hbW8POzg61atWCnZ0d+Hw+cnJykJubCxMTExgYGJS7doFAgFevXuHVq1cAPnbdWllZQSgUQiAQ4M2bN0hPT+dqyc7ORm5uLvLy8sDn86GlpQVdXV1YWFjAysoKNWvWhJ6eHnR1dVFQUICkpCQkJSUhJSUFAoEAAoEAPB4P1atXh5mZGbZv3y7T710lgmXr1q04c+YMduzYwU1buHAh3rx5U+a4iJg4WM6cOQNbW9vPLrtgwQJERETAwcEBhw4dAgD4+PggLy8PJ06cQP369XHlyhX07t0bWlpaOHHiBJydnUusJzk5GbGxscjOzsaHDx/QsmVL2NnZybDninPz5k1069YN1tbWuHHjxmf/ANXR+/fvce7cOZw6dQpJSUkIDAxE7969MX/+fBw8eBC1a9dGaGgoJk2ahPz8fInn2tvb47fffkNAQABGjBiBY8eOQUtLC0VFRbC0tMSRI0dgaWmpoD2rXAKBAHv37sXu3bvx/PlzqZ6jp6cHd3d3NG7cGP7+/vD29oauri6AjwH17NkzXLhwAQ8ePIC5uTlq1aqFWrVqwdLSEhYWFtDW1kZiYiLi4+Px5s0b5Ofno6CgAE+ePMHVq1eRm5tb5ra1tbW5H/F7XBwo/31dK8rQ0BA1a9aEtrY2eDweNDQ0oKurC319fWhpaSEvLw95eXnIzs7Gu3fvkJmZKdfty6I8vUKfUon+jQ8fPnBvNDFdXV25X5pj3LhxuHDhAh4+fIjBgwfD09MTeXl5CAwMRP369QEAvr6+6N+/P3bt2oXJkycjKipK4g25evVqhIeHS3zbNTY2xs6dO+Hl5SXXeuVtyZIluHv3Lpo2bYqHDx8C+Dho/7WFCvDxNevSpUuJgxYWL16MuLg4PHnyBD/++CMAoH///vj5558RFRWFHTt24NmzZxgyZAgcHR2RmJiIatWqYd++fZg5cyZu3LiB4OBg/P7779DT0wNjDLGxsYiJiUF6ejoCAwMRGBhY4v1elfLy8nDnzh3cvHkTt2/fRmFhIfdt38TEhPv2nZCQgEePHiE+Ph5CoRA8Hg8FBQXcTdUsLCzg6+uLpk2bok6dOsjIyMDbt28hEAiQmpqK1NRUxMXF4dWrV7h9+zZu376NtWvXQldXF6amptDU1ERhYSHS0tIqtD9OTk6wtrbG+/fv8f79e+Tk5HA/hYWFn73qta6uLqpVqwYjIyNUq1YNNWvWRM2aNVGjRg2YmpqiWrVqKCoq4loVycnJ+Pfff5GSkgLGGIyMjKCnp4fMzExum9LS0NCAubk5rK2tYWVlBQDc701LSwvm5uYwNzeHmZkZ97oYGhrC0NAQ+vr6EIlEKCoqQm5uLtLS0pCamgqBQID8/Hzk5+dDU1MTtWvXRp06dWBnZ8ftF4/Hw7t375CRkSHz71ylWywCgQCLFy/+7HPL02IBPo6pBAUFSZww92nXGPDxSr+tW7dGWloaZs+ejY4dOyIzMxO//fYbd993X19fmJmZIS0tDTdu3ICBgQF27NgBHx+f8u5+lbh48SK+//77EtPFLTXyP8+fP0eHDh2QlZWF8ePH4+eff+a6aIuKirBlyxYsWbIEOTk50NTUxO7du+Hr64t3796hY8eOSE5O/uz6jY2N0bVrV0yaNAnVq1evlH3IysrCgwcPkJWVBU9PT9jY2CA9PR2bN2/G5s2bZT4Kks/nIzAwEP369UPr1q2l6uJKT09HdHQ0rly5gvPnz+Px48cS801MTNCyZUv4+PggKysLL168QHJyMtLS0pCWlob8/HzUrl0b9erVg62tLXR1daGtrQ1ra2u0aNEC1tbWpW5XJBJxwVJYWMidJM3j8aCnpwc9PT2Zv1QJhUJoaGhw7wvGGDIzMyEQCFBcXAzGGEQiEfLy8vDhwwcUFhZCT08P+vr6MDIyQo0aNWBiYqKyX+pUIlguXbpUYvCzvGMs0gYLAMTFxaFr167IyspCy5YtsWfPnhLLHDlypNSr+daoUQPLly/nzocoLi7G+PHjERUVBT09PezYsQPNmjWTqo7KkpOTg5cvX3LdeAUFBQgMDERSUhL69u0LPp+Pq1evwt3dHWvXruX+OMj/JCcn4/Xr12WOPb158wabNm2Cj48PAgICuOkJCQlYsGAB94EoEong6OgId3d36Ojo4ODBg1xrsUaNGli4cGGZtyMQCoW4cuUKTp8+DXt7e/Tp0weGhoYllsvIyMCDBw8QExODmJgYPHz4EElJSRLL2NjYICMjAx8+fADw8XYDzZs3h4+PD0xMTPD+/XtkZWUhMzMTGRkZyM7ORq1atdCwYUO4uLhwrS8+nw8dHR2Zfqdi4m/3QqEQAGBra/vZgBKJRCr7AayuVCJYCgsL0bZtW/Tv3x9DhgzB5cuXERoaisOHD39x7EKWYAGAW7duISIiAtOnTy91HIUxhqlTp+Lo0aPQ0dGBrq4uGjRogN9++w01a9aUWFYoFOLnn3/G3r17YWFhgbNnz8LExETqWuQlMTERW7duxd69e5GTk4N+/fphzpw52LBhA8LCwuDo6IjTp09DW1u7ymsj//P06VPMmDED165dA/DxqLzRo0ejUaNGAD6G0759+/D333/j9evX3POMjIzQt29fWFlZITMzE2/evMHdu3e5++V8SkdHB+7u7qhWrRru3LmD9+/fAwBatWqFcePGKW3LmqgGlQgW4GMr4tdff8XTp09hYWGByZMnS3WWtKzBIm9CoRDdu3fH7du30atXLyxfvrxKt79p0ybMmjWLe8zn8yEUCuHu7o7ExETk5+fjzz//RIsWLaq0LlI6kUiErVu3Yt68edwgcqNGjaCtrY1bt25xy9nb26Nz5864ffs2rl+/Xuq6dHR00LBhQzRo0ADu7u5wc3ODi4sL9wVCJBIhNjYWWlpadOUJIhcqEyyyUpZgAT5+02zfvj3y8/Oxbds2BAYGcvNEIhGSkpKQnp6OJk2ayLX76datW+jRoweEQiH69u2L4OBgMMYQEhKCFy9eAACCgoKwZs0auW2TyEdycjK2bNmCv/76iztKyMDAAJ07d0bv3r3h7e3NvVcePHiAyMhIMMZgYmICU1NTNGjQAA0aNKhw9xQh5UHBUsXWrl2LOXPmwNLSEqNHj0ZiYiJiY2Px6NEj7ogR8VFD8giXjIwMtGvXDq9evcLIkSMxc+ZMbl5WVhamT5+O+Ph4bNu2jTvyhCgf8SHvjDG0a9eu3OdEEFKVKFiqmFAoRNeuXXH37t0S86ysrPDu3TsUFhZi0KBBmDdvHgDg7t27SE9PR8OGDct1/gNjDMHBwTh58iQ8PDwQGRlJ4yeEkEqnEuexqBM+n4+IiAjMnTsXpqamcHR0RL169VC/fn3UrFkTZ86cwQ8//IDt27cjNjYWSUlJEpcLsbS0RMeOHfHLL798sXtjx44dOHnyJIyNjbFmzRoKFUJIlaAWixI6f/48hg0bxg3a2tnZwc7ODg8fPkR2djYAwNvbG5s2bSrzPIfU1FS0bt0a2dnZWLNmDYKCgqqsfkLI141aLEqoVatW2LdvH65fvw5/f3+4ubmBx+NBJBLh7t27GDFiBG7evImgoCAsXrwYjo6O3BmzYjNnzkR2djbat2+Pzp07K3BvCCFfG2qxqKDU1FQMHjwYMTEx3DRdXV34+/tj1KhRePv2LYYPHw4DAwOcP3++zDOPCSGkMlCLRQVZWVkhMjIS8+bNw507d/Dy5UtkZmbixIkTOHHiBHdl3WnTplGoEEKqHAWLijI0NMSCBQu4x2/evMGWLVuwbds2ZGVlwcPDA4MGDVJghYSQrxV1hamZnJwcnDt3Dr6+vpV2AUNCCPkcarGoGUNDQxqsJ4QoFF0SlBBCiFxRsBBCCJErChZCCCFyRcFCCCFErihYCCGEyBUFCyGEELlS+8ONxffN/vQWroQQQqRjaWkJTc3yRYXaB4v4kvP9+/dXcCWEEKJ6ZDm5XO3PvM/Pz8ejR49gbm4OPp+v6HIIIUSlyNJiUftgIYQQUrVo8J4QQohcUbAQQgiRKwoWQgghckXBQgghRK4oWAghhMgVBQshhBC5omAhhBAiVxQsSmjv3r1wd3eHh4cH9xMVFaXosuQmOjoa3t7e3OPCwkLMnj0bPj4+8PHxwYIFC7hL8aiy/+5nfn5+idc1ODhYgRXKLjo6Gv369UPjxo3RqlUrrFy5EowxtXsty9pPdXotAeDcuXPo3LkzPDw80LZtW/z5558AKvC3yYjSmTlzJlu6dKmiy6gUR44cYV5eXuybb77hpi1evJj169ePZWRksNTUVNalSxe2bt06BVZZcaXt571791iLFi0UWJV85OTksKZNm7IdO3aw4uJi9uzZMxYQEMB27typVq/l5/ZTXV5Lxhh7+fIlq1+/Prt06RJjjLEnT56wBg0asAcPHsj8elKLRQnFxMTA1dVV0WXI3bJly7Bx40aMHj1aYnpUVBRCQkJgYmICS0tLjBw5Evv371dQlRVX1n6qy+uampoKT09PDBgwAHw+Hw4ODmjbti3u3r2rVq/l5/ZTXV5LALC2tsa1a9fg5+cHkUiEjIwM8Pl8GBoayvx6UrAomaKiIsTFxSEyMhJ+fn5o27Yt1q9fD6YGV97p168fIiMj4e7uzk17//49BAIBHB0duWl16tTBixcvUFhYqIgyK6y0/QQ+BotAIEDnzp3RvHlzjBs3DmlpaQqqUnZ169bFqlWruMeFhYW4ePEiXFxc1Oq1/Nx+qstrKWZoaIicnBzUr18fQ4YMwcCBA1GjRg2ZX08KFiWTnp6Ohg0bonv37jh79iz++OMP7N69Gzt37lR0aRVmYWFRYtqHDx8AAHp6etw0PT09rh9bFZW2nwCgr68PLy8vbNu2DceOHYOOjg5GjRpVxdXJV2FhIUJDQ6GtrY0OHToAUK/XUuzT/ezfv79avpZ6enq4f/8+9u3bh7/++ov7zJHl9VT7y+arGgsLC+zatYt77ObmhoEDB+LkyZMYOHCgAiurHOI37adv1Ly8PAAfP4jVyYwZMyQeT5s2Dc2aNUNKSkq5L0uuDAQCAX788UcAwJYtW6Ch8fF7qrq9lv/dT319fbV7LQGAz+eDz+ejYcOG6NmzJ6KjowHI9npSi0XJxMbGIiIiQmJaQUEBtLW1FVRR5apWrRrMzc2RlJTETUtKSoKDg0O5L9WtzBhjWLZsmcR+irsTdHR0FFWWzOLi4tCjRw84ODhg+/btMDU1VcvXsrT9VLfX8tq1a+jdu7fEtKKiIhgbG8v8elKwKBkDAwOsX78e+/fvh0gkQnR0NHbu3ImePXsqurRKExQUhFWrVuHdu3dIS0vD2rVr0bVrV0WXJVc8Hg8xMTFYtGgRsrOzkZWVhblz58Lf3x/m5uaKLq9cMjIyEBwcjE6dOiEsLEziS486vZZl7ac6vZYA4Orqin///Rfbt2+HUCjE7du3ERUVhZ49e8r8etL9WJTQ5cuXsXTpUjx79gympqYYNmyYWt0B88aNGxg5ciTu3bsH4GOLbOHChThx4gSKi4vRpUsXTJkyReVvzPbf/Xz37h3mzp2Lq1evQigUwt/fH7NmzUK1atUUXGn5bNmyBWFhYdDT0wOPx+Omt2jRAuHh4WrzWn5uP2fPnq0Wr6XYw4cPMW/ePMTHx8Pa2ho//vgj2rVrJ/PfJgULIYQQuaKuMEIIIXJFwUIIIUSuKFgIIYTIFQULIYQQuaJgIYQQIlcULIQQQuSKgoWovICAADg7O5f687W7ceMG97v4/fffkZKSAmdnZ8TFxZVYNjIyEj4+Pl9c57Jly7h1njt3rjLKJipONa+zQMh/hIaGonv37oouQ2mdPn0apqamyMzMrPC6QkJCMGDAAPj5+VW8MKKWKFiIWjAwMFDJy2lUlerVq8PAwEAuwWJgYAADA4OKF0XUFnWFEbW3cuVKDB8+HMOGDYOXlxd3m+dNmzahdevW8PDwwPfff4/79+9zzykuLsbixYvRtGlT+Pj4YP369ejevTsiIyMBAAMHDsTChQsltvPfrqHPrX/q1Kn49ddfMXXqVHh4eKBZs2YS9/4QiURYtWoVWrVqBQ8PDwwaNAiJiYl48+YN3NzccO3aNW5ZxhgCAgLw999/y+13FhkZWWrXojpeYZvIHwUL+SpcunQJ3t7e2Lt3L/z9/fHnn39i586d+PXXXxEVFQV/f38MHjwYKSkpAD6GUVRUFBYtWoStW7fi7NmziI+Pl3p7X1o/APz999+wtrbGwYMHMXToUKxYsQKPHz8GAERERGDHjh2YPn06oqKiYG5ujpCQEJiZmaFp06Y4evQot5579+7h7du3+Pbbb+X02wI6dOiAy5cvcz/Lly+HpqYmhg8fLrdtEPVFwULUQlhYGDw8PCR+bt++zc3X09NDSEgIHB0dUb16daxbtw6hoaHw9/eHg4MDRo4cCS8vL+zatQuMMezduxejR49Gy5Yt4erqiiVLlkAoFEpdz+fWL1arVi2MGzcOtWrV4m7/+vDhQzDGsGfPHowePRrt2rWDg4MDZs+ejbZt2yIrKwtBQUE4efIkioqKAAD//PMPWrduDUNDQ6nr69WrV4nf1+zZs7n5urq6MDc3h7m5OYqLizFnzhyMGTMG/v7+Um+DfL1ojIWohREjRiAoKEhi2qd3crSxseGuUJubm4tXr15hxowZmDVrFrdMYWEhtLW1kZ6ejvT0dLi5uUk839LSUqpavrR+sVq1akk8z8DAAMXFxcjIyEB6ejoaNGjAzTM2NsbUqVMBAG3btsWvv/6Kq1evws/PD8ePH8fvv/8uVW1iK1asgIODg8S0kydPYuPGjRLTCgoKMGbMGHzzzTcqf4dEUnUoWIhaMDU1hb29fZnzP70Bk7jlERYWJhEewMdv6rq6ugA+jnN86nM3WysuLpZ6/Z9bH2MMWlpa3P9LY2BggDZt2uD48ePQ1NTkLtteHlZWViV+X2ZmZiWWmzlzJj58+IBFixZJXDqekM+hrjDy1RHfGS8tLQ329vbcz7Zt23Dp0iUYGBjAxsYGDx484J6TkZGBly9fco+1tbWRm5vLPU5OTpZ6/V9iZGQEMzMzbrwF+Hh7WF9fX+52sUFBQTh37hzOnDmD9u3bc2EkT1u3bsXp06exatWqcnWzEUItFvJVGj58OFavXo2aNWuifv36OHToEP766y/s2LEDwMdzNZYuXQo7OzvUrl0bixcv5m4/CwD169fHvn370KVLFxgYGJS4k+KX1v8lQ4YMwZo1a2Brawt7e3usXr0ahoaGcHFxAQD4+fmBz+fj77//xubNm+X4m/no2rVrWLx4MebPnw8TExMIBAJuHh3WTb6EgoV8lQYNGoT8/HyEh4fj7du3qF27NlasWAFPT08AQN++fZGZmYlff/0VBQUFGDZsGHcnSAAIDg5GUlISgoODYWpqinHjxkm0aL60/i8JDg5Gbm4uZsyYgdzcXHh6emLdunVcePH5fHz33Xc4e/YsvLy85Pib+ejgwYMoKirCpEmTSsyLjY2V+/aIeqE7SBIiJR8fH0yZMkVpzvAfN24cateujQkTJpS5zI0bNzBo0CDcvXtX7ic1Ojs7Y+3atWjdurVc10tUH7VYCFExN2/exOPHj3HhwoVSWxSlSU9PB2NMLmMlubm5+PDhQ4XXQ9QXDd4TomKOHDmCFStWYNKkSbCzs5PqOYGBgVi6dKlctr9+/Xq6Thj5LOoKI4QQIlfUYiGEECJXFCyEEELkioKFEEKIXFGwEEIIkSsKFkIIIXJFwUIIIUSu/g8NB6/Dx4aIoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the oscillatory component on a linear-linear scale\n",
    "plt.plot(freqs, psd_osc[1, :], 'k', lw=2)\n",
    "plt.xlim(1, 30)\n",
    "sns.despine()\n",
    "plt.title('Oscillatory component, chan = ' + chan[1])\n",
    "plt.xlabel('Frequency [Hz]')\n",
    "plt.ylabel('PSD log($uV^2$/Hz)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the oscillatory + aperiodic component on a linear-log scale\n",
    "# psd_combined = psd_aperiodic[1, :] + psd_osc[1, :]\n",
    "# plt.plot(freqs, psd_combined, 'k', lw=2)\n",
    "# plt.fill_between(freqs, psd_combined, cmap='Spectral')\n",
    "# plt.xlim(1, 30)\n",
    "# plt.yscale('log')\n",
    "# sns.despine()\n",
    "# plt.title('Aperiodic + oscillatory component, chan = ' + chan[1])\n",
    "# plt.xlabel('Frequency [Hz]')\n",
    "# plt.ylabel('PSD log($uV^2$/Hz)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the fractal component (1/f)\n",
    "\n",
    "If `return_fit=True` (default), `yasa.irasa` [fits an exponential function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) to the aperiodic power spectrum and return the fit parameters (intercept, slope), the $R^2$ of the fit, and the standard deviation of the oscillatory component.\n",
    "\n",
    "The aperiodic signal, $L$, is modeled using an exponential function in semilog-power space (linear frequencies and log PSD) as:\n",
    "\n",
    "$L = a + log(F^b)$\n",
    "\n",
    "where $a$ is the intercept, $b$ is the slope, and $F$ the vector of input frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02-Jun-22 14:35:07 | WARNING | The evaluated frequency range ends after the resampled Nyquist frequency (52.63Hz). Decrease the upper band (30.00Hz) or decrease the maximum value of the hset (1.90).\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>Chan</th>\n",
       "      <th>Intercept</th>\n",
       "      <th>Slope</th>\n",
       "      <th>R^2</th>\n",
       "      <th>std(osc)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>F4-A1</td>\n",
       "      <td>3.351799</td>\n",
       "      <td>-1.297857</td>\n",
       "      <td>0.961306</td>\n",
       "      <td>1.012285</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CZ-A2</td>\n",
       "      <td>3.272062</td>\n",
       "      <td>-1.192108</td>\n",
       "      <td>0.855980</td>\n",
       "      <td>4.525508</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Chan  Intercept     Slope       R^2  std(osc)\n",
       "0  F4-A1   3.351799 -1.297857  0.961306  1.012285\n",
       "1  CZ-A2   3.272062 -1.192108  0.855980  4.525508"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "freqs, psd_aperiodic, psd_osc, fit_params = yasa.irasa(data, sf, ch_names=chan)\n",
    "fit_params"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}