notebooks/11_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 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/manual/io.html#importing-data-from-eeg-devices\n",
    "\n",
    "Note that during loading, [MNE automatically scale the EEG data from microVolts (uV) to Volts](https://martinos.org/mne/dev/manual/io.html#id3). However, YASA requires the data to be in microVolts (uV), and therefore we need to reverse this scaling by multiplying our data by 1e6 (1 V = 1,000,000 uV)."
   ]
  },
  {
   "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() * 1e6  # Convert data from V to 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 [10_bandpower.ipynb](10_bandpower.ipynb)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(401,) (2, 401)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEgCAYAAAD2c3e8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3gU5fr/8fem914hkZBEIKFID6FKFykCPxSQIsWjFCmKX0ERsCGiSAnCCWKJVGmJqOABBUR6h1ATQhFCCSUJkL7Znd8fOe5xDSXZbDJJuF/XxXXOPjM7c0/i5rMz88zzaBRFURBCCCFUYqF2AUIIIR5vEkRCCCFUJUEkhBBCVRJEQgghVCVBJIQQQlUSREIIIVQlQSREGcnKymLu3Ll07tyZevXq0a5dOz755BPu3r0LwKBBg6hZs+Z9/02aNOmR29+/fz81a9Zk6tSp912+bds2evXqRYMGDejSpQtr16416/EJYSqNPEckROnLyMigf//+ODg4MHr0aKpXr87FixeZOXMmdnZ2LFu2jJycHLRardH7Nm/ezAcffMC8efN45plnHrqPd999lwMHDnDr1i127dqFnZ2dYdnx48fp378/b775Ju3atePw4cNMmTKF2bNn07Fjx1I5ZiGKSs6IhCgDs2bNQq/X891339G6dWsCAwNp1aoVixcv5syZM6xbtw43Nze8vb0N/+7du8esWbPo06fPI0MoLy+PTZs2MWLECLRaLf/5z3+MlsfFxdGsWTOGDBnCE088Qc+ePXnuuef44YcfSvOwhSgSK7ULEKKyy8vL46effuLNN980OksB8Pf3Z8mSJVSvXr3Qe9588018fX2ZPHnyI/exbds2MjIyaNOmDS1atGDdunX07NnTsLxfv37o9Xqj92g0GnJyckpwZEKYh5wRCVHKLl++TEZGBnXr1r3v8oYNG+Lu7m7UNmfOHBITE5k9ezYODg6P3MePP/5Iw4YN8fDwoGPHjhw4cIDLly8blteoUYNatWoZXt+4cYMNGzbQuHFjE49KCPORIBKilN25cwcAZ2fnIq2/e/duvv32WyZMmEB4eHiRtr99+3bDvZ527dphaWnJunXr7rt+RkYGo0ePxtPTk0GDBhXxKIQoPRJEQpSyv852/gqkh0lNTeWtt96iRYsWDBkyxGjZjz/+SIMGDQz/oqOjAfjll1/QarV06tQJADc3N5o2bUpcXFyhy3FpaWkMHTqUa9eusXjxYpycnMxwhEKUjNwjEqKUVatWDTc3N44fP069evUKLZ8xYwbe3t68/PLLTJ48Gb1ez8yZM9FoNEbrtWvXjqeeesrw2tXVFSgIKIAOHToYlun1ehRFYefOnbRu3RqAlJQUhg4dSnZ2NsuWLSMoKMjchyqESSSIhChlFhYWPPfccyxbtow+ffpga2trWJacnMyqVat44403WL58Odu2bWPRokV4eXkV2o6Tk1OhM5grV65w+PBhxowZYzgjAsjPz2fgwIGsW7eO1q1bk5GRwbBhw9BqtSxfvpwqVaqU3gELUUwSREKUgVGjRrF9+3YGDx7MmDFjqFatGgkJCXz22WfUrFmTpk2b0rdvX3r37k14eDg3b940er+1tTVubm6Ftvvjjz9ia2vL4MGDcXFxMVrWq1cvVq1aRVpaGtHR0Vy+fJmYmBisra0N27eysirUUUKIsiYPtApRRm7fvs3ChQvZtm0bt27dwsfHh06dOjFy5Eh+/fVX3n777Qe+t2nTpixdurRQ+7PPPstTTz3FjBkzCi07f/48zz77LO+88w5fffUVKSkphdapXbs2sbGxJTswIUpIgkgIIYSqKuWluTt37jBhwgRyc3OxtbVl1qxZ972sIYQQQn2Vsvv28uXL6dy5M0uXLqVz586sWLFC7ZKEEEI8QKU8Ixo0aBA2NjZAQe8hK6tKeZhCCFEpVMq/0H89wZ6YmMjKlSvve5NXCCFE+VApL80BHD16lAkTJjB37lzDg39FkZ+fT3JyMvn5+aVYnRBCiL9UyiBKSkpi6tSpREdHExwcXKz3Xr9+nfbt23P9+vVSqk4IIcTfVcpLcwsWLCA7O9swq2W7du0YOnSoylUJIYS4n0oZRHPmzFG7BCGEEEVUYS7NxcfHExkZadR25swZ+vbtS/369enevTvx8fEqVSeEEMJU5T6IFEVhzZo1hgEb/5KXl8eoUaPo0qULBw4cYMSIEQwfPpyMjAwVqxVCCFFc5T6IoqKiWLlyJSNHjjRq379/P1qtliFDhmBtbU3Xrl0JDQ1l48aNKlUqhBDCFOU+iPr160dsbCx16tQxak9KSiIkJMSoLTg4mMTExLIsTwghRAmV+yDy9fW9b3tWVhZ2dnZGbfb29mRnZ5dFWUIIIcyk3AfRgzg4OJCbm2vUlp2djYODg0oVCSGEMEWFDaKQkBAuXLhg1Hb+/HlCQ0NVqkgIIYQpKmwQRUREoCgKMTExaLVaNmzYQEJCAh07dlS7NCGEEMVQYYPIxsaGxYsXs2nTJpo2bUp0dDQLFizAw8ND7dKEEEIUg8zQ+g/Jycm0b9+eLVu2EBAQoHY5QghR6VXYMyIhhBCVgwSREEIIVUkQCSGEUFWlHH1bPB7WrVvHsWPH8PLywsfHhxYtWhAYGKh2WUKIYpIgEhXSpUuXGDduHH/va+Pn58euXbsKjbghhCjf5NKcqJBWrlyJoii4+T2Jf3gnXNx9uX79Ot9//73apQkhikmCSFQ4+fn5rFq1CgDv8G741+uJZ3h3oGB23ry8PDXLE0IUkwSRqHC2bNlCSkoKbl5VcfIuGNLJLaA+bt4BXL16lTVr1qhcoRCiOCSIRIWzbNkyAJwCm6HRaADQaCxwCy0Y3umLL74wmkRRCFG+SRCJCuXKlSts27YNSytrPIIijJa5BzbC1bMKly5dIi4uTqUKhRDFJUEkKpTvv/8eRVHwDmqIla2T0TKNhQVuoR0AmD17tsxNJUQFIUEkKgytVsuKFSsAcA5sdt91PKo1wd3nCS5fvsyiRYvKsjwhhIkkiESFsXHjRq5fv17QScGnxn3X0VhY4lW7F1Bwr+jq1atlWaIQwgQSRKLC+OabbwBwCWpl6KRwP86+NfENbkx2djbTp08vq/KEECaSIHoAmR2jfDl27BgHDx7E1t4R92pNH7m+V/hzWFvb8MMPP7B///4yqFAIYSoJogfYsWOH2iWIv/n6668B8ApujqX1o4fwsXXyxCfsf925hRDllwTRA3zzzTdyVlRO3Lhxgx9//LGgV1z1VkV+n0dwazQaDTt27ODevXulWKEQoiQkiB7gzJkzbN26Ve0yBAUPsGq1WnyC6mPr5FXk91nbOeMdUJO8vDy2bNlSihUKIUpCgugh5syZI2dFKrt9+zaLFy8GwKVa0c+G/mLjXQco6HEnhCifJIgewN7RmSNHjrB9+3a1S3mszZw5k7t37+JTrd4Du2w/jFvV+gBs3bpVHnAVopyq9EEUExNjuNFdHG7VmgMFT+jrdDpzlyWKID4+nhUrVmBhaYlX7Z4P7bL9IDaOHnj5h5Cdnc0ff/xRClUKIUqqUgfRtGnTWL58uUnvdavWBHtHFw4dOsTYsWNlEM0ypigKU6ZMQVEU/Gu1w87Fz+Rt2frI5TkhyrNKHUStWrVi5MiRJr3X0sqWgIhh2Nja88MPPzBs2DC5tFOG4uLiOHjwII7ObnjVeqZE23L97+W5X3/9Vb5QCFEOVeog6tChQ4ne7+xTg+ptxmLv6MLWrVsZNGgQOTk5ZqpOPIhOp2P27NkAeIV1w9LavkTbs3Pxxd0nkDt37rBnzx5zlCiEMKNKHUTm4OhRjaDW43B08WDPnj28+eab0pOulP3yyy9cuHABF3ffQlM9mMrBty4Aixcvlt+fEOWMBFER2Lv6Exg5Ahtbe+Li4pgzZ47aJVVaiqKwcOFCANyCn0ZjYWmW7XoEt8LWzoGtW7eyefNms2xTCGEeEkRF5OAeQGDEEDQWFnz++efExsaqXVKltGvXLo4dO4aDkyse1SPNtl1re1d8a3cFYOrUqXK/T4hypNIHUe/evRk+fLhZtuVapS6BDfsA8NZbb3Hu3DmzbFf8z4IFCwDwCGmDhZWNWbftGdIad59qJCcnM3/+fLNuWwhhugoVRPHx8URGGn9LPnPmDH379qV+/fp0796d+Pj4Uq3BK/Rp/EIjyc7Olm7dZnb8+HH++OMPbGzt8Agu/igKj6KxsMS7bsEXiX//+98kJSWZfR9CiOKrEEGkKApr1qxh2LBhRn/48/LyGDVqFF26dOHAgQOMGDGC4cOHk5GRUar1+Nbrg7ObN0ePHmXevHmluq/HyaxZswDwDmmJla1jqezDyTsEvxotyMvLY/z48fJFQohyoEIEUVRUFCtXriz0TND+/fvRarUMGTIEa2trunbtSmhoaKk/uGhpY49fwwFoNBqioqI4dOhQqe7vcbB7925+++03bGzt8ahRsm73j+JTpxdOrl4cOXKEuXPnluq+hBCPViGCqF+/fsTGxlKnTh2j9qSkJEJCQozagoODSUxMLPWanH1qUKV2R3Q6HSNHjuT27dulvs/KSq/X8+GHHwLgU6sj1nYupbo/KxsH/BsNNHyRkInzhFBXhQgiX1/f+7ZnZWVhZ2c8SZq9vX2Z9YjyDuuKh38oV65c4ZVXXpHLPCZav3498fHxOLl44PVkuzLZZ8EXic7o9XrGjBnDnTt3ymS/QojCKkQQPYiDgwO5ublGbdnZ2Tg4OJTJ/i0sranSdBiOzh7s3buXadOmlcl+K5OcnBxmzJgBgGetZ83eU+5hvMOexcOvOsnJyYwZMwa9Xl9m+xZC/E+FDqKQkBAuXLhg1Hb+/HlCQ0PLrAYbezcCmg3Hysqa7777jnfffZdr166V2f4rutmzZ3PlyhXcfZ7AI6hZme7bwtIK/8ZDsXNwYsuWLYbOEkKIslWhgygiIgJFUYiJiUGr1bJhwwYSEhLo2LFjmdbh6FmdgEb9Afj222+JjIzkjTfeYP369Zw5c4a8vLwyraei+O2331iwYAEaCwu86/RBY1H2/znaOnkR0HQoGgsL5s2bx4YNG8q8BiEedxU6iGxsbFi8eDGbNm2iadOmREdHs2DBAjw8PMq8Fo/qzajVaRJ+IU3I1+lYtWoVo0aNon379tStW5dffvmlzGsqz5KTkxk7diwAAfV64ORTdmex/+TiF0ZA/V4AjB8/nvPnz6tWixCPI40iI0AaSU5Opn379lRvNwlrB9MCLefeDdIvHUDJuk5W2hXupF7HxtaW71euJCLCPIN4VgQ6nY6ff/6ZrKwsmjVrRlBQEFqtlmPHjvHee+9x9OhRfILqUzXiZTQadb8TKYrC1QPfkXJ+P40aNSIuLg5LS/OMcyeEeDgrtQuojOycffD777hmiqKQcmwNV8/8zpAhQ4mLi6VWrVoqV1j6jh07xsSJEzl+/Lihzc/Pj/T0dMNUGk5u3vg1GKB6CAFoNBp8679Axq0kDh06xOLFixkxYoTaZQnxWFD/L0Alp9Fo8K3XB5/qDbl79w4DBgzg7NmzapdVahRFYcaMGXTt2pXjx4/j5OqFf0hj7B1duH79Ojk5Obj7BFKl1tNUa/FaqY2gYAorGwf86vcDYOann5bJ82hCCDkjKhMaCwuqNBqMLi+T61cS6NatG/Pnz6dTp05ql2Z2X331FV988QUWFpZUCW+PV61nsLS2w0/Rk3vvBpY2jljbOatd5gO5VqmDf42WXEvcyeuvv86PP/4ol+iEKGVyRlRGLKxsCIgcgV9IEzIyMhg6dCizZ8+uVM+u7NixwzBCQlDkUHzr9sTSuuCBY43GAjsXv3IdQn/xqVswBNDRo0dLfbgoIUQJgkiv15OUlMS+ffs4ePAgFy9elJkvH8HSyhb/xkN4omEvNBoNn3/+OS+//DL37t1Tu7QSu3jxIiNGjECn0xFQtwtugQ3VLslkltb2eNYoeARg3rx5lerLghDlUbEvzR04cIAlS5awe/duMjMzDe0ajQYXFxdatWrFiy++SMOGFfcPUWnSaDR41eiIjXNVkg/EsGnTJrp168Znn31mGCniiSeewNvbW+1SH2n//v0sX76cY8eOkZSUhKIo+FR7Cu/wrmqXVmLuQc24nbiZ06dP89tvv1XKy6hClBdF7r598eJFpk6dyrVr1+jQoQOtW7cmJCQENzc3FEUhNTWVM2fOcODAATZu3EhgYCAffPAB1atXL+1jMCtzdN8uqpx7N7i6/2vSb142are2tub5559n1KhR5fLnp9PpmD9/Pp9//rnhbMHC0hLvwNr4NXwJSxt7lSs0j5tnt3H50BqeeuopNmzYgEajUbskISqlIgfR888/z+jRo3n66acfua6iKPz2228sWrSItWvXlrTGMlWWQQSg0+Zw8+QPZN8+h4WlNRqNBbevX0BRFCwsLBgzZgxvvfVWqddRVCkpKYwfP54//vgDjUaDf3hHnKvUx96tKhaW1mqXZ1b6/DwS//MeWRnpLFu2jLZt26pdkhCVUpGDSFGUYn8jNOU9aivrILqfnLsppJ3bSsrZ3ej1OlavXk2LFi0My/Pz88nNzSU/Px9ra+tSG+Q1MzMTa2trbGxsyM7OJjo6moULF5KVlYWDkytVGg/GxS+sVPZdXtxM+JXLR+IICAigW7duPPXUU7Ru3Ro3Nze1SxOi0ihyZ4W/B8o777xDQkJCoXVSU1Np3779fd8jis7OxRf/Bv2pUudZACZNmmR4CDQuLo7atWtTo0YNwsPDqV27Nt99951Z968oClFRUYSHhxMcHEyjRo2IiIhg1qxZZGVl4RvUgOD2Eyt9CAF4BLfCxd2X5ORkoqOjGTlyJG3btuXIkSNqlyZEpWFSr7nY2Fj69evHzz//bNSu1+u5evWqWQoT4FWzA25eVTl//jxffPEF33//PWPGjCEjIwNrG1ts7R3Jy8tj6tSpHDhwACgIkc8//5xnn32W9evXF7sno06nY/LkycycOZP8/HxAw/Xr17l9+zbuvkE82XY8VZv9Cxv7x+OMwNLajurt3iKk1Uiq1O2Kh18wN27c4P/9v//H+vXrC62fn5/P+fPn//uzE0IUhUljzdWqVYu3336buXPn8sILLzBx4kQsLCy4desWrVq14vTp06VRa5koD5fm/u7ejbOc3ToHKysrwx+3wPo98a5V0Isr5XgcV07+ir+/Pxs2bODjjz82ui/XuHFj3n//ferXr//IfWVlZTFu3Dg2btyIlbUNTzQdjGvVp8jLSkenzcberUq5GI5HTYpex/Vja7iW8AcAXbp0oUWLFtSuXZutW7eyevVqUlJSqFatGq+88gp9+/bF3r5ydN4QorSYHES7du3ixo0bjB49moCAAObOnYuiKLRs2VKCyMyuHV7BtcSdADzR6Hm8nvzfTXNFr+PSjihuXzuLk5PTf8+W7PAMbcPdS/vIykgHoE+fPkyaNAkvLy9++uknli1bhqOjI2PHjqVJkyacOnWKkSNHkpSUhJ29I4GRr+Ds86Qqx1veKYrC7aTfuXx4HYpS+BkjaxtbtHkFEza6u7vTpUsXOnbsSKtWrSSUhLgPk4IoLCyMnTt34unpSWpqKuPHj+fSpUvMmDGDYcOGSRCZWX5eFrdP/4itewhuTzQptDwvM5VzWz8lO/Mu9o4uPNF8BI6eQei02dxO/JXrp7eQn6/F3t4eNze3QhP3RUZGcujwYfJyc3HzCsC/yVDsXf3L6vAqrNyMW9xLOYM2/Ty5d65g6xaIU0AETl7BpF+J5+75rdy+ds6wvourK998/TWRkZEqVi1E+VOiMyJPT0+g4L7CjBkzWL16NVqtVoJIBZm3LnAveR9uwW2xc/E1WpabcYvbp3/k+rmDALh5VcWleht0OencTNxKXm5BRwj/mq3xrdu7TKfrrswURSE7/Qr3rh0n+8ZxUq9fxNHRkVWrVtGgQQO1yxOi3DApiL744guGDx9e6DJDXFwcsbGxLF261GwFlrWKGkRFkZV2GZ02GyfvUMO9Hm3OPW6f24mdaxXcAp5SucLKS9HruXZ4GdeT9uLq6so333xDSkoKu3fvxsPDgzFjxpRaN3whyrsiB1FSUlKRNxoaqt5smyVVmYNIqEvR67iy/xtuXCzc9btmzZp8+eWXhs9ORXwGTwhTFXmsuW7duqHRaAwfkL9/UP6eZRqNpkJfmhOitGgsLKnSZAgoelKvnsazSg0sXauTeeUwCQkJPPvss3Tp0oWzZ8+SkJBAZGQkX375pZwpiUqvyGdEV65cMXqtKArdu3fnyy+/pEqVKkbLqlatar4Ky5icEYnS9tdH7q8vcjptDtePfk/Kuf2F1o2MjOS7777D0bH8TCAohLkV+YzoQeHi5+dXoYNHiLL2z0tultZ2VGn8Eo7+9cnLSisYt8/KluS9i9mzZw8DBw5k6dKlODk5qVSxEKWrUj6dmJ+fz4QJE3jxxRcZN24cubm5apckxENpNBrcAurjU6Mtzj41cPSoRrVWY3Fy9WT//v306NGDkydPql2mEKWiUgbR5s2b8fPzY8WKFdSrV48ff/xR7ZKEKDY7Zx+qtRqHi6c/CQkFU8xHR0fLRH2i0qmUQXT06FFatmwJQIsWLdi/v/C1dyEqAlsnL6q3fYsqtdqQl5fHhx9+SMuWLfn3v/9Namqq2uUJYRZFvke0fPnyQm06nY7169fj7u5u1D5gwICSV1YCGRkZhpu7jo6ORjPJClHRWFrZ4le/L/Y+4dyIX8uff/7JRx99xMyZM6lduzb16tUjIiKC7t27Y2lpqXa5QhRbkYPo66+/LtTm5eVFbGysUZtGo1E9iJycnAzhk5mZKT2ORKXgWqUuLn61uXvtJBmXd3Pjz3iOHj3K0aNHWbJkCUeOHOH9999Xu0whiq3IQbR169bSrMOs6tWrx65du4iMjGTXrl1FGnlaiIpAY2GBa9W6uFati2/DLLLTLpN5+wIpJzfw1VdfERERwbPPPqt2mUIUS5GDqFOnTrRo0YIWLVrQrFmzct2V9JlnnuH333+nX79+uLu7M2fOHLVLEsLsrGwccPatibNvTaysbbl0aA1vvPEG4eHhBAUFGa176dIlTp48iaOjI87OztSqVUtGAhflRpEfaL1w4QK7d+9m586dHDp0iODgYFq2bEnz5s2pX78+FhaVo9+DPNAqKiJFUQqGD7pwiNq1a7Ns2TJ8fHyAgoksJ0x4k7y8/z3G4OPjy5Il31G3bl21Sn6gJUuWsGrVKgYOHMjzzz+PlVWRvy+LCsqkQU91Oh1Hjhxh9+7d7Nq1i4sXL9K4cWPDGVO1atXMXmh8fDyvvvoqe/bsMbSdOXOGadOmkZCQQGBgINOnT6devXol2o8EkaiodHnZnN/2KffSUnBycmLs2LHcu3eP+fPnA+AdUAtLS0uy7t3ibmoKDg4OREdH0759+wduc/369XzzzTf069ePvn37lvoXzpUrV/Lmm28aXoeEhDBp0iS6dOkiY+9VYiYF0T/du3ePPXv2sGvXLvbs2cPmzZvNURtQ8E1v7dq1zJw5E4CDBwumMsjLy+OZZ55h8ODBDBgwgM2bN/Pee++xbdu2El02lCASFVlu5m1uHl/HjYtHDW0WFpYENuyDZ2gbAPQ6LdePrOR60l4sLCwYPXo0/fv3N/oCmZ+fz/Tp0/nyyy8NbREREcycOZMnnzSeMHH9+vVs3LgRFxcXXF1dqVGjBt26dcPBwQG9Xs+aNWuIjo7G2tqayMhImjdvTsuWLQt1ItqwYQMjRoxAr9fj82Rrcm4ncjf1OgDPPfccM2bMwNXV1ew/M6E+swRRaZo3bx7bt2+na9eu/Pvf/zYE0c6dO3n77bfZsWOHYd3+/fvTq1cvXnjhBZP3J0EkKoO7105x+/R6sjPSCWjyEi7+4UbLFUXh5qkNJB/faGhr0qQJTz75JDY2Npw8eZIDBw5gaWmJV2hr7l05QlZGOtbW1kRFRdGjRw8AfvrpJ0aMGFFo/87OzvTs2ZMjR45w4sSJQsvt7e3p0qUL3bp149atWxw7dow1a9aQl5dHwFPd8Qnrgl6XT+r5XVw/sZ683BwCAgKYP38+TZs2NfNPS6itSEF0/vx5UlJSqFu3rtHZxrZt22jbtu1D3llyKSkp+Pr6sm/fPkaPHm0IopiYGH7//XdiYmIM606ePBl7e3veffddk/cnQSQqC0VRQNGjsXjws0X3UhLJSN7LrT+Pos3LMVrm4OxGQNNhOHmHkp+Xxc0TP3AtcScWFhbMnz+f0NBQnnvuOXJycvAL64CNow/5uRnkpZ7m1pWzhu04uXriFdYNa3s3Mm8loU09y80rifetp0p4B3zr9jK6DJdz7wbXDy0h9fp5rKysmDlzJv369SvhT0eUJ4+8C7h8+XKWLl1K9erVOX36NJMnT6Zjx44AREVFlXoQ+fr63rc9KysLOzs7ozZ7e3uys7NLtR4hKgqNRgOahz/g6uxbA2ffGvjUe4G710+jy81Er9ei0VjgFlAfa/uCS2FWNg74NeiPpZ0ryfEbGDNmDO7u7gUhVKMF/vX+Hh7P4JWeTEbyYTQ2TngGtzTM+uvsWxMAn4xbpF86gDYtEUs7dywcq+LkFYyDZ1Che0F2zj5Ua/06did+5OqpX5kwYQJ//vknb731ltw3qiQeGUTff/89sbGxODg4cPnyZcaNG0dycjJDhw5Fzat6Dg4OhQYzzc7OlrlbhDCBpbUd7oEPn75co9HgE94VNBqSj/3M7du38fR/Er+nXigUCA5uATi4BTxwW7ZOXviGdwG6FKk+jYUlfvV6Ye3kzeWDq4iKiuLs2bO88847BAcHF2kbovx6ZBDpdDrDH/fAwECWLFnCuHHjuH79uqpBFBISYnRZDgouIfbs2VOdgoR4TPiEPYuVjSO5qefwqvP/sLC0LrN9ewa3xMrOncv7vuGXX35h06ZN9OrVi6CgIA4ePEh8fDzVq1fnxRdfpEePHjKqSgXxyCDy8vLi9OnThIWFAQXD5yxatIh33nmHxMT7X+ctCxERESiKQkxMjKHXXEJCguGyoRCi9HiEtIGQNqrs27VKbew6vk1q4q+knNvNunXrjJanpaVx+PBhpk2bxtixYxk1atQju50fOnSIbdu24awSEAQAACAASURBVObmhre3Nw0aNOCJJ54ozcMQf/PIzgrXr1/H0tISb2/vQssOHTpEo0aNSq24v/tnZwWAxMREpk2bxpkzZwgICOCdd94hMjKyRPuRzgpCVBy5GbdJv7gDvU6HnUd1HNyfIOPGWbKu7udmcgIAbdq0ISoqCi8vr/tuY9myZbzzzjvodDpDm42NDatXr6ZJkyZlchyPu2J1346Li6NXr16lWY/qJIiEqBzuXD3J1UNLyc68i4+PL/PmzaV169aG5Xq9nunTpxMdHQ2AT0gkVjZ25N29xq0rZ/D29mbjxo1UqVIFRVE4ePAg/v7+BAQ8+N7XgyiKQlZWFmlpaTg4OODhIX9b/q5Yj0lPmTKFAQMGqHpJTgghisK1Sm1C2k/Cs2pNbtxIoX///kyZMoXMzEw2bNhgmGjQ0tKSoIhBBDQZhN9TzxPYYjReAeHcvHmTl19+mcOHD9OvXz969uzJ008/zapVq4p8fzw7O5uJEycSGhpKjRo1iIiIoEGDBkyfPp2srKxS/glUHMU6I0pMTOTDDz/kyJEjDBw4kDFjxlS6m4FyRiRE5aLoddxK+JWrJzag0+mwtbU19Li1d3SlapOXcPGrZfSe/NwMLmybxb30G4Y2K2sb8rV5APTq1Ytp06YZ3bK4c+cOSUlJ+Pn54e/vz9mzZxkxYoThi7uVtQ32ji5k3LmNoihUqVKFjz76iM6dO5f2j6DcM2lkhQ0bNvDZZ5+h1+uZNGlSpRp2XoJIiMopK+0y1w8vI/3mZVzcfXALbotH9UjDM06F1k9P5vy22eTna/Gr+TReNTtz9+oJrhxZbXj4NzAwkLCwMC5evGh0pcjOzg69Xk9eXh6unlXxb/KSoTt75u2L3Dy+htTrF4CCUPvwww8LTTD6ODF5iJ/s7GwWLlxITEwMjRs3ZurUqVSvXt3c9ZU5CSIhKi+9Tkt2ejIO7k88dMSJv+RlpQNg4+BmaMu5m8LtU3GkXUskL/d/o1FYWlnj5lmF7Mw7ZGUUvM+vRgt86/XB0srWaLuKXs/tc9u5fvwn8vJy8PHxYfjw4Tg7O2Nvb0/z5s1NuhdVUZV4rLmLFy8yffp09u7dy7Bhwxg5cmShEQ8qEgkiIURRKHo9OXevk33nCjaOnji4BxqeqcrPy0Kfn2cUYPeTc+8GKUdXcvtKglG7vb09H3/8cYnGzfyn7OxsvvnmG65evYqXlxeenp6kpqZy7tw5Ll26RLVq1WjRogUtW7akatWqj9zewYMHmT9/Pmlpabz00ks899xzJk/ZYbZBT7du3cr06dNRFKVCzeb6TxJEQoiypOj1pP65D+29q6DPR5t1m5uXCgaK7dOnDx9//PEj78Xn5+cTHx/PzZs3ad68Oc7OzkbLd+/ezZtvvsmff/5ZpJrq1atH7969ee655wzzWgHcvHmT/fv3s2LFCn7//Xej91SvXp2dO3cWafv/ZHIQXb58mYSEBBITEw3/++eff6IoCqdPnzapmPJAgkgIoSZFUUi7uPe/96JyCQkJITo6mvBw4xHUU1NT+eWXX9i8eTN79+4lIyMDAFtbW9q3b0+DBg24du0a586dY/v27QC4+zyBQ5Um6LUZaPKzUCztsHb0wcbRk5y7V8lPP0/a1TPk5vyvR5+jo6PhGay/B5mNrT3eNZ7Gyt6DtKQt3E29zpUrV0w65mIF0dSpU0lMTCQxMdEwuGhAQABhYWHUrFmTWrVqUatWrSKd1pVXEkRCiPIg+85Vrh38jvSbl7GxtWXa1KkEBQURHx/P/v372bFjB/n5+Yb13byqYG3rdN+RzS0trfCv0wWvGh2xsHz45TN9fh53rp4g+9ohbl0+QX6+1rDMxtYOd79QrN2C8QhuiZVtwWwMil5H6oXd7P5hlknHWqwgeuGFFwxhU6tWLWrWrCndt4UQopTo8/O4cTyWqwl/FFpmYWGJd2A4tj71cPELw8ax4O9VXlYa6ZcPo89JR2Prio2jJ44e1QzLi0NRFHTabPJzM1B0+di5+D60k8dPnz9X7H1AEcaa+7vVq1ebtBMhhBDFZ2Flg1+Dfth6hnLn3G9Y2Thg5VQVa5cAXPzDsbZzLvQeGwd3fGo+ePr34tBoNAX7tCndWQ1M6uLwww8/ULt2bYKDg7G0fHQXSCGEEKZzf6Ix7k80VruMUmNSEH322WekpqZibW3Nk08+SVhYmOFfrVq1ZE4gIYQQRWZSEO3atYuUlBROnjzJqVOnOHbsGOvXryc/Px+NRsOpU6fMXacQQohKyrSnjyiYwtvX15d27doBkJKSwhtvvEGPHj3MVpwQQojKr1ijbz+Mr68v06dPL/SQkxBCCPEwJgVRamrqfduDgoLkspwQQohiMenSXPPmzfH19SUsLIzw8HDCw8MJDg4mPj4erVb76A0IIYQQ/2VSEK1evZpTp05x6tQptm/fzldffUVeXh6WlpZMnDjR3DUKIYSoxEwKonr16lGvXj3Da51Ox5UrV3B1dcXV1dVsxQkhhKj8inyP6O7duw9cZmlpyRNPPFEohO7cuWN6ZWYUExPD119/rXYZQggh7qPIQTRw4EC+/PJL7t2798h1U1NTWbhwIQMHDixRceYwbdo0li9frnYZQgghHqDIl+ZWrlzJ3Llzad26NQ0bNqRVq1Y8+eSTuLu7oygKqampnDlzhgMHDnDw4EF69erFypUrS7P2ImnVqhVPPfUUaWlpapcihBDiPoocRI6OjkyePJlXXnmFVatW8fPPP3P69Gl0Ol3BhqysCA8Pp02bNnz44Yf4+vqWWtH/tHbtWlasWGF47ePjQ3R0NAAdOnQgNja2zGoRQghRPMXurODt7c1rr73Ga6+9hl6vJz09HY1Gg7u7e2nUVyR9+vShT58+qu1fCCGE6Uwe4gfAwsICDw+Zs0cIIYTpTAqiWrVqodFo7rvM2toaX19funbtypgxY2SaCCGEEA9lUhB98MEHzJs3jzFjxlC/fn0Ajh8/TlRUFC+88AJPPvkkX3zxBYqi8Prrr5u1YFP07t1b7RKEEEI8gEljzX399ddMnz6dfv36GaYNf/755/n444/ZsGEDzzzzDO+//z4//PCDyYXFx8cTGRlp1HbmzBn69u1L/fr16d69O/Hx8SZvXwghRPlgUhDduHGDwMDAQu1+fn5cu3bN8P9NeaBVURTWrFnDsGHDjMaty8vLY9SoUXTp0oUDBw4wYsQIhg8fTkZGhimHIIQQopwwKYgaNmzIZ599Rnp6uqEtLS2Nzz//nIYNGwKwefNmgoKCir3tqKgoVq5cyciRI43a9+/fj1arZciQIVhbW9O1a1dCQ0PZuHGjKYcghBCinDDpHtFHH33Ev/71L1q3bk3VqlXR6/Vcu3aNkJAQoqKi+P3335kzZw7z588v9rb79evHuHHj2Ldvn1F7UlISISEhRm3BwcEkJiaacghCCCHKCZOCyN/fn/Xr17N3714SEhKwsrLiySefNNzTcXR05Pfffzepa/eDHoTNysrCzs7OqM3e3p7s7OziH4AQQohyw+TniCwtLbGwsMDKygq9Xo9WqyU/Px8rK6tSebbIwcGB3Nxco7bs7GwcHBzMvi8hhBBlx6QgunnzJqNHj+bUqVNUrVoVRVG4evUq1atXJyYmBk9PT3PXSUhICDExMUZt58+fp2fPnmbflxBCiLJjUmeF6dOnY2lpyZYtW9i0aRObN29my5YtuLm5MWPGDHPXCEBERASKohATE4NWq2XDhg0kJCTQsWPHUtmfEEKIsmFSEO3YsYN3333X6H6Or68vEydO5I8//jBbcX9nY2PD4sWL2bRpE02bNiU6OpoFCxbIEENCCFHBmXRpztbW9r5D/Gg0GsNo3CUVERHBwYMHjdpq1KhRLqaWEEIIYT4mnRG1aNGCGTNmcOvWLUPbrVu3+OSTT2jZsqXZihNCCFH5mXRG9NZbb/HSSy/Rtm1bqlSpAsDVq1epWbMmkydPNmuBQgghKjeTgsjb25sff/yRHTt2kJSUhJ2dHSEhITRv3tzc9QkhhKjkihxESUlJhdoCAwONxpz7a53Q0FAzlCaEEOJxUOQg6tatGxqNBkVRDG3/7LCgKAoajYbTp0+br0IhhBCVWpGDaMuWLaVZhxBCiMdUkYOoatWqpVmHEEKIx5RJ3beFEEIIc5EgEkIIoSoJIiGEEKoq9nNEqamp/P7775w9e5bMzEycnJwICwujTZs2uLi4lEaNQgghKrFiBdG3337LvHnz0Gg0BAQE4OzsTEZGBkuXLsXKyorXX3+dwYMHl1atQgghKqEiB9G6deuYN28eb731Fr179zaaLTU3N5fY2Fg+/fRTfH196dy5c6kUK4QQovIpchAtXbqUN998kxdffLHQMltbW/r3709WVhZLliyRIBJCCFFkRe6scPHixUeOrN2uXTvOnTtX4qKEEEI8PoocRDk5OTg5OT10HWdnZ+7cuVPiooQQQjw+itV9+36T4QkhhBAlUaxec7GxsTg4ODxweWZmZokLEkII8XgpchBVqVKlSNN0+/v7l6ggIYQQj5ciB9HWrVtLsw4hhBCPKZNmaK0I7ty5w4QJE8jNzcXW1pZZs2bh5uamdllCCCH+oVidFTIyMvjhhx/IyMgACibCW7x4MYMHD2b8+PHEx8eXSpGmWL58OZ07d2bp0qV07tyZFStWqF2SEEKI+yhyEF26dIkuXbrwwQcfkJaWBsCMGTOYPXs2Xl5euLi48NJLL3H48OFSK7Y4Bg0aRI8ePQDIz8/HyqrSnvwJIUSFVuS/zvPmzaN27drMnTsXOzs70tLSWLFiBc888wyzZ88GoFq1anzxxRd88803pVbwP61du9bobMfHx4fo6GicnZ0BSExMZOXKlSxdurTMahJCCFF0RQ6iPXv2sGjRIsMYczt37kSn09GzZ0/DOi1atGDhwoXmr/Ih+vTpQ58+fe677OjRo0yZMoV58+bh6upapnUJIYQomiIH0b179/D09DS83rt3L1ZWVjRt2tTQ5ujoiF6vN2+FJkpKSmLq1KlER0fLNOdCCFGOFTmIqlatSlJSElWqVEGn07F9+3aaNGmCvb29YZ29e/cSEBBQKoUW14IFC8jOzmbSpElAwTh4Q4cOVbkqIYQQ/1TkIOrduzcfffQRr732Gnv37uXWrVtMmzbNsPzgwYPMnTuXQYMGlUqhxTVnzhy1SxBCCFEERe41N3z4cDp06MD06dPZtm0bb7zxBh07dgTgww8/ZODAgTRq1Ijhw4eXuKj4+HgiIyON2s6cOUPfvn2pX78+3bt3L1ddxYUQQphOoyiKUtKNJCQkoNfrCQsLK9F2FEVh7dq1zJw5Eyg4ywLIy8vjmWeeYfDgwQwYMIDNmzfz3nvvsW3btkeOCF5cycnJtG/fnurtJmHt4GHWbQshRGX20+fPmfS+Ip8RZWdnM2XKFJo2bUrLli157733DA+21qxZs8QhBBAVFcXKlSsZOXKkUfv+/fvRarUMGTIEa2trunbtSmhoKBs3bizxPoUQQqiryPeIoqKi2LZtGy+//DIWFhYsX76c9PR05s6da7Zi+vXrx7hx49i3b59Re1JSEiEhIUZtwcHBJCYmmm3fQggh1FHkINq0aROzZs2iWbNmADRt2pQBAwag1WqxtrY2SzG+vr73bc/KyjI8v/QXe3t7srOzzbJfIYQQ6inypbmUlBSCg4MNr+vWrYuiKNy+fbtUCvs7BwcHcnNzjdqys7MfOjeSEEKIiqHIQaTT6bC0tDS81mg0WFtbo9VqS6WwvwsJCeHChQtGbefPnyc0NLTU9y2EEKJ0FWv0bbVERESgKAoxMTFotVo2bNhAQkKCofu4EEKIiqtEU4XrdDrWr1+Pu7u70XoDBgwwT3X/ZWNjw+LFi5k2bRrz5s0jICCABQsW4OEh3auFEKKiK/JzRO3atSvaBjUatmzZUqKi1CTPEQkhhGlMfY5IpgoXQgihqgpxj0gIIUTlJUEkhBBCVRJEQgghVCVBJIQQQlUSREIIIVQlQSSEEEJVEkRCCCFUJUEkhBBCVRJEQgghVCVBJIQQQlUSREIIIVQlQSSEEEJVEkRCCCFUJUEkhBBCVRJEQgghVCVBJIQQQlUSREIIIVQlQSSEEEJVEkRCCCFUZaV2AaUlIyOD8ePHk5mZSXh4OO+++y4ajUbtsoQQQvxDpT0jWrNmDW3btmXlypVkZWVx7NgxtUsSQghxH5X2jGjo0KHodDry8/NJTU3FyclJ7ZKEEELcR4UPorVr17JixQrDax8fH6KjowHQ6/V069YNV1dX/Pz81CpRCCHEQ2gURVHULqK0LV++nBs3bvD6668/ct3k5GTat29P9XaTsHbwKIPqhBCicvjp8+dMel+lvUe0aNEitm7dCoC9vT0WFpX2UIUQokKr8JfmHqRnz55MnDiRb7/9FicnJ2bOnKl2SUIIIe6j0gaRr68vMTExapchhBDiEcrl9ar4+HgiIyON2s6cOUPfvn2pX78+3bt3Jz4+XqXqhBBCmFO5CiJFUVizZg3Dhg1Dq9Ua2vPy8hg1ahRdunThwIEDjBgxguHDh5ORkaFitUIIIcyhXAVRVFQUK1euZOTIkUbt+/fvR6vVMmTIEKytrenatSuhoaFs3LhRpUqFEEKYS7kKon79+hEbG0udOnWM2pOSkggJCTFqCw4OJjExsSzLE0IIUQrKVRD5+vretz0rKws7OzujNnt7e7Kzs8uiLCGEEKWoXAXRgzg4OJCbm2vUlp2djYODg0oVCSGEMJcKEUQhISFcuHDBqO38+fOEhoaqVJEQQghzqRBBFBERgaIoxMTEoNVq2bBhAwkJCXTs2FHt0oQQQpRQhQgiGxsbFi9ezKZNm2jatCnR0dEsWLAADw8ZC04IISq6cjmyQkREBAcPHjRqq1GjBitXrlSpIiGEEKWlQpwRCSGEqLwkiIQQQqhKgkgIIYSqJIiEEEKoSoJICCGEqiSIhBBCqEqCSAghhKokiIQQQqhKgkgIIYSqJIiEEEKoSoJICCGEqiSIhBBCqEqCSAghhKokiIQQQqiqXE4DoSadTgeANjtd5UqEEKJiSU5Oxs/PDyur4kWLRlEUpZRqqpAOHjzIgAED1C5DCCEqpC1bthAQEFCs90gQ/UNOTg4nTpzA29sbS0tLtcsRQogKRc6IhBBCVDjSWUEIIYSqJIiEEEKoSoJICCGEqiSIhBBCqEqCSAghhKokiIQQQqhKgkgIIYSqJIgqgbVr11K7dm0aNGhg+BcXF6d2WWYTHx9PZGSk4XVeXh5TpkyhadOmNGvWjEWLFqlYnfnc7zjr1Klj9HsdNmyYihWWzK5du+jduzcNGzakY8eOfP/99wDcvXuXsWPH0qhRI1q1asW6detUrtR0DzrGy5cvU6tWLaPf5eTJk1Wu1nTbtm2je/fuNGjQgA4dOhiO0+TPpiIqvPfff1/57LPP1C7D7PR6vbJ69WqlUaNGSqNGjQzts2bNUgYMGKCkp6crly9fVjp37qzExcWpWGnJPOg44+PjlebNm6tYmflcvXpVadCggbJ582ZFp9Mpx44dU5o0aaL88ccfyuuvv66MGzdOycrKUk6dOqVERkYq+/btU7vkYnvYMf7yyy/K888/r3aJZpGSkqLUqVNH+f333xVFUZQTJ04odevWVU6cOGHyZ1POiCqBkydPEhYWpnYZZhcVFcXKlSsZOXKkUXtcXBwjRozA1dWVgIAAhg8fbvhGVhE96DhPnjxJrVq1VKrKvK5cuUK3bt3o2LEjFhYW1KtXj6ZNm7Jnzx42bdrEuHHjsLe3JywsjOeff55Vq1apXXKxPegYDx8+XKl+lz4+PuzZs4c2bdqg1+tJT0/H0tISR0dHkz+bEkQVnE6nIyEhgfXr19OyZUs6duzIl19+iVIJRm7q168fsbGx1KlTx9B29+5dbt68SWhoqKGtevXqJCYmqlGiWdzvOAFOnTpFamoq3bt3p3nz5owdO5aUlBSVqiyZxo0b88EHHxhep6enc/DgQXx8fNBoNAQFBRmWBQcHV8jf54OOMTw8nFOnTnH27Fk6d+5My5YtmTx5Mnfv3lWx2pJxcnIiOzubunXrMmzYMAYMGICHh4fJn00JogouNTWVOnXq0LNnT7Zu3Wr4dr1ixQq1SysxX1/fQm1ZWVkA2NnZGdrs7e3Jyckps7rM7X7HCQXH1bBhQ7777jv+85//YGdnx+jRo8u4OvO7d+8eI0eO5KmnnqJ27drY2Nig0WgMyyv67xOMj7F9+/a4uLjQokUL1q5dS1xcHNeuXWPKlClql1kitra2HDlyhLVr17Ju3TqWLFkCmPbZlPmIKjhvb2+WLVtmeB0WFsbAgQPZvHlzpZzOwt7eHoDc3FxDW3Z2Ng4ODmqVVGrefvtto9eTJk0iMjKSa9eu4e/vr1JVJXPhwgVGjRpFaGgos2bN4ty5c+Tl5aEoiiGMKvrv85/HaGFhwZw5cwzLnZ2def3113nxxRfJz88v9kjV5YWFhQU2NjbUrVuXF154gRMnTgCmfTbljKiCO3v2LFFRUUZtWq0WW1tblSoqXa6urnh7e3P+/HlD24ULF4wuB1QW8+bN49y5c4bXWq0WoML+bg8cOMALL7xAhw4diIqKwtbWlmrVqqEoCsnJyYb1zp8/X2F/n/c7xuzsbD799FNu3bplWE+r1WJlZVUhp5rZv38/vXv3NmrLy8vDxcXF5M+mBFEF5+Liwrfffsvq1avR6/WcOHGCpUuXFvoPpTLp0aMHCxYsIDU1leTkZL7++mt69Oihdllml5CQwCeffMLdu3e5e/cu06dP5+mnn8bDw0Pt0ort0qVLvPrqq4wdO5YJEyYYzn4cHR3p0KEDn3/+OZmZmZw5c4Y1a9ZUyN/ng47R3t6e3bt3M3v2bHJycrhx4waff/45vXr1MrokWVGEhYWRkpLCt99+i06n4/Dhw6xbt44+ffqY/NmU+YgqgT179vDZZ59x4cIF3N3dGT58eKW6LLdv3z5Gjx7NwYMHgYJT/08++YRNmzah1+vp27cv48ePr5Af6r/753GmpaXx0UcfsXPnTnQ6HW3atGHq1Km4urqqXGnxzZgxg5iYmEKXaV588UVeeeUV3n//fXbt2oWNjQ3/+te/GDx4sEqVmu5hx9i/f38++ugjDh06hIWFBV27dmXixIkV9uz25MmTfPTRRyQmJuLv78/YsWPp1KmTyZ9NCSIhhBCqkktzQgghVCVBJIQQQlUSREIIIVQlQSSEEEJVEkRCCCFUJUEkhBBCVRJE4rHTrl07atased9/+fn5apenqn379hl+Fh988AHJycnUrFnzvgNXxsbGEhER8chtzpo1y7DNbdu2lUbZooKrmIMcCVFCb7zxxn1Hn6io436Z22+//Ya7uzvp6ekl3tbIkSN56aWXaNmypRkqE5WRfOrEY8nR0RFvb2+1yyi3PDw8cHR0NEsQOTo64ujoaIaqRGUll+aE+If58+fz8ssvM3z4cBo1amSYdv3rr7+mbdu2NGjQgP79+3P06FHDe/Lz85k5cybNmjUjIiKCL7/8kt69exMbGwvAoEGDmDlzptF+/nmp6mHbnzRpEu+99x5vv/02DRo0oF27dixcuNCwXKfTsWDBAp5++mkaNGjAoEGDOHv2LDdu3CA8PJw9e/YY7btjx45mnZI7Njb2vpc6Bw0aZLZ9iMpLgkiI+9ixYweNGzdm9erVtGnThu+//54lS5Ywbdo04uLiaNOmDS+99BKXL18GCsJr/fr1zJw5k5iYGLZu3crZs2eLvL9HbR9g7dq1+Pr6GgaYnDdvnmHo/QULFrB06VLeeecd4uLi8Pf359VXX8XT05NmzZrx888/G7Zz9OhRUlJS6NSpk5l+WvDss8+yc+dOw7+5c+diZWXFyy+/bLZ9iMpLgkg8lj755BMaNGhg9G/fvn2G5XZ2drz66quEhITg4eHBokWLmDBhAk8//TRBQUGMGDGCRo0asWLFChRFYfXq1YwcOZI2bdoQFhbGrFmz0Ol0Ra7nYdv/S2BgIOPHjyc4OJhRo0bh5ubGyZMnURTFMNV4p06dCAoKYurUqXTq1Ik7d+7Qo0cPfv31V8M0Ej///DNPP/00zs7ORa7v+eefL/TzmjZtmtHPy9vbG29vb/Lz8/nwww8ZPXo0bdq0KfI+xONL7hGJx9Krr75aaHj6v8+UGhAQgIVFwfe0zMxMrl69ypQpU4z++Obl5WFjY0NqaiqpqanUrl3b6P1+fn5FquVR2/9LtWrVjN7n6OhIfn4+aWlppKamUrduXcMyJycnJk2aBBRchnvvvffYuXMnrVu35pdffuH9998vUm1/iYqKMprOG2Dz5s189dVXRm25ubm89tpr1K9fn5EjRxZrH+LxJUEkHkvu7u6F/rD/3d+H59fr9UDBWVR4eLjRenZ2dkZTI//d30Pkn/7eTfxR23/Y9hRFwdraGuCBQ+07OjrSvn17Nm7ciI2NDXl5ebRu3fqBtd2Pv79/oZ+Xp6dnofWmTJlCZmYmn376aYWflkOUHbk0J8QjODs74+3tTUpKCtWqVTP8W7ZsGTt27MDR0ZGqVasadS5IS0vjypUrhtc2Njbcu3fP8Prv934etf2i1Ofp6cmpU6cMbTk5OTRv3pxjx44BBZMJ/vHHH/z2228888wzDw1JU8XExPDbb7+xYMECnJyczL59UXnJGZEQRfDyyy+zcOFCfHx8qFu3Lj/99BPLly9n6dKlALzyyivMnj2bwMBAqlevzqxZs8jLyzO8v06dOixZsoSePXvi5OTEJ598YhQGj9r+owwZMoSFCxdStWpVgoKCiI6OxtnZmbCwMABatmyJlZUV69atK3Q5zRz27NnDrFmz+Pjjj3Fzc+PmzZuGZdJNXjyKBJEQRTB48GBycnL47LPPuHXrFkFBQURFRdGoUSMA+vXrR3p6Ou+99x65ubkMHz6cI0eOGN4/bNgwzp8/OoJuWAAAARZJREFUz/Dhw3F3d2fs2LFGZ0yP2v6jDB8+nMzMTN59910yMzNp2LAhixYtMoSdpaUlXbp04ddff6VJkyZm/MkUWL9+PVqtlv/7v/8rtCwhIcHs+xOVi8zQKkQpiYiIYOLEifcdwUENY8eOpVq1akyYMOGB6+zbt4/Bgwdz+PBhsz+EWrNmTaKjo2nbtq1ZtysqPjkjEqKS279/P6dOnWL79u1GzxM9TGpqKoqimOVeT2ZmJllZWSXejqi8pLOCEJXchg0biIqK4v/+7/8IDAws0ns6dOjA7NmzzbL/f//73zLOnHgouTQnhBBCVXJGJIQQQlUSREIIIVQlQSSEEEJVEkRCCCFUJUEkhBBCVRJEQgghVPX/ASkP+9pu6+HUAAAAAElFTkSuQmCC\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": [],
   "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": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEgCAYAAABvi64JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde1xM+f8H8Nd0kVKpsKjWNV2kKEREQqGotChbWuSS+7ouWnate1jkskiyKpalUi6bRa5FIoVQqVQkkaikppnP7w9f52e2UJlpuryfj8c8dudzzpzzOjNm3p3b58NjjDEQQgghNUBG2gEIIYQ0HFR0CCGE1BgqOoQQQmoMFR1CCCE1hooOIYSQGkNFhxBCSI2holOPxcTEQE9PD8uXL5f4urZt2wYnJ6dqv379+vUYN24cAOD69evQ09NDUVGRuOI1SHw+H0FBQWJZVnBwMHr16iWWZUnS27dvsWXLFgwZMgTGxsYYOHAg1q1bhzdv3nDzjBs3Dnp6ehU+Fi9e/MV1fOl7FRkZiZEjR8LExATDhg3D0aNHxbZ99QEVnXosLCwM7dq1w8mTJ/Hu3TuJrmvixInw8/MTy7JMTExw5coVKCkpiWV5DdWJEyfg4+Mj7Rg1prCwEM7OzoiOjoaXlxdOnjyJFStW4MqVK5g4cSL3Hdi2bRuuXLki8vhQQAYMGPDF9Xzue3Xnzh3MmjULDg4OOH78OKZOnYoVK1bg33//Ffv21lVUdOqp0tJSREREwNPTE3w+H//8849E19ekSROoq6uLZVmNGjVCixYtwOPxxLK8hqqh3fe9ceNGCIVC/Pnnn+jfvz++/fZb9OvXD76+vnjw4AGOHTsGAFBTU0OLFi24R0FBATZu3IhRo0Zh6NChn13Hl75XISEh6N27N8aPH482bdrA0dERDg4OCA0Nldh21zVUdOqpyMhIFBYWwtLSEn379uW+cMD7QyUODg7YvXs3zMzM0Lt3b/z+++8QCoXcPPHx8XBxcYGRkRFsbGzg6+vLTQ8ODsbIkSOxYMECmJqaYufOneUOr925cwfu7u4wMTGBhYUFNmzYAD6fz02/dOkSRowYAWNjY8yYMQMFBQXctP8eXnvy5AmmT58OU1NT9OnTB6tWrRJZ1sfy8/OxePFimJmZwczMDAsXLuSWXVpaCh8fHwwcOBBGRkYYO3Ys4uPjudeOGzcOf/zxBzw9PWFsbAwbGxtcv34dhw4dQr9+/dCzZ0+sWrWKm3/x4sVYvnw55s+fj65du8La2honTpwQyRMeHs5t55AhQxASEsJN27ZtG2bOnIn169fDzMwMFhYWWLVqFQQCATdPaGgohgwZgq5du2LkyJG4cOFCpV5//fp1LFmyBPn5+dDT08P169crfL/+6+DBg9z6Ro0ahdjYWJHpfn5+sLCwQLdu3TB79myRzy0sLAwjRoxAly5dYGpqCk9PT+Tm5nKfaa9evRAaGoqBAweie/fu8PT0xIsXLyrMsW3btk8eAqtoW0pLSxEeHg43Nzc0btxYZFrr1q1x4MAB2NraVvi6BQsWoGXLlvDy8vri+/O57xUAuLi4YMGCBSJtPB5P4kca6hRG6qXp06ez77//njHG2LFjx5ienh7LyMjgnhsaGjJnZ2eWmJjIzp8/z8zMzNi2bdsYY4y9ePGCde/enW3fvp2lpaWxCxcusAEDBrDdu3dzr9fV1WXLly9n6enp7MmTJ8zHx4eNHDmSMcZYamoq69q1K1uxYgVLSUlh58+fZ3379mXr1q3jphsaGjIfHx/26NEjtnv3bqanp8fc3NwYY4xdu3aN6erqssLCQlZSUsJsbGzYxIkTWWJiIouNjWVWVlbs999/r3C7XV1d2ciRI9mtW7dYYmIic3R0ZPPnz2eMMebl5cX69+/PLly4wFJSUtiyZcuYiYkJy8nJYYwx5ubmxoyNjdmRI0dYeno68/DwYD169GAeHh4sKSmJex+joqIYY4z99NNPzNDQkC1fvpylpKQwf39/pq+vz65du8YYY+z48ePM0NCQBQUFsbS0NBYQEMAMDQ1ZZGQkY4wxHx8fZmhoyJYuXcoePXrEDh8+zPT19dk///zDGGPs0qVLrHv37uzEiRPs8ePH7NChQ8zIyIjdunXri68vKSlh+/fvZ2ZmZuz58+espKTki/9mjh49yoyNjdnff//N0tPTmbe3NzM1NWV5eXncZz516lSWlJTErl+/zszMzNjGjRsZY4zdvHmTGRoaspCQEJaVlcWioqKYpaUl++2337jP1MDAgDk7O7M7d+6w6Oho1rt3b276fxUWFrLnz59X+KhoW1JSUpiuri67c+fOF7fzY+vWrWOGhobs3r17lZr/c9+riuTk5DATExO2c+fOKuWqz6jo1EP5+fnM0NCQ+fv7M8YYe/XqFevcuTPbvHkzY+z9l0VfX59lZmZyr/H392d9+vRhQqGQbd26lY0fP15kmWFhYczMzIx7va6uLnv+/Dk3/eOis27dOjZ8+HAmFAq56SdPnmSGhoasqKiIeXt7MycnJ5HlT5w4scKiExkZyQwNDdnLly+5eS9dusSCgoLKbXdSUhLT1dVl9+/f59oSEhLYjh072OvXr5mBgQE7ffo0N00gEDBbW1uugLm5ubHJkyeLZNbV1WXp6elcm5WVFTtw4ABj7H3RsbGxYQKBgJs+bdo0NmfOHMYYYyNHjmQrVqwQybh8+XLm7OzMvWempqYiP6KOjo5s06ZNjLH3BfS/P1Y///wzmzVrVqVef+zYMe4zqwwnJye2evVq7rlQKGTe3t4sPT2d+4H9+HNYunQpmzRpEmOMsXv37rHg4GCR5f3222/M3d2dMfb/n2lcXBw3fc2aNdx78bVu3rxZ7rP6kqtXrzI9PT22b9++Ss3/pe/VfxUUFLBRo0axwYMHs4KCgkrnqu/kpL2nRcTv9OnT4PP5sLGxAfD+GLaZmRlCQkIwe/ZsAO8POWhra3OvMTY2xosXL/Dq1SukpKQgJiYGJiYm3HShUIh3797h1atXAABFRUW0aNGiwvUnJyeja9euIudkunfvDj6fj8ePHyM5ORmGhoYirzE2Ni53KAcAUlJSoKWlBQ0NDa6tX79+Fa43JSUF8vLy0NPT49qMjIxgZGSE+Ph4CAQCkW2SkZGBiYkJkpOTubY2bdpw/6+oqAgA+Pbbb7m2Ro0aobS0lHtuYmICGZn/P0ptbGyM8PBwLs/48eNFMnbv3h0nT57knmtqaqJRo0bcc2VlZe7QYXJyMuLj47Fnzx5uOp/PR/v27Sv1+qp69OiRSF4ej4eFCxcCAG7evAkVFRWRz0FVVRUZGRkAgM6dO6Nx48bYvn07UlNT8ejRIyQnJ6N79+4i62jXrl2lsu7atQu7d++ucJqvry969Ogh0vbhfOLr168rta15eXlYtGgR+vbtW+4zCgsLwy+//MI9nzp1Kjw9Pb/4vfr438GrV68wZcoUZGdnIzAwEMrKypXK1RBQ0amHwsLCAACDBw/m2oRCIRhjuHLlCgBAVlZW5DUfztfIyMigrKwMNjY2+PHHH8stW0VFBQCgoKDwyfX/95g68P8ntT+sh/3nJLe8vHyFy/pU+6fm/dTFB5/KyxgTOZdV0fo+/jH5r4rexw/zf+p9+NL6PhAIBJg/fz6srKxE2uXk/v9rW5X350s+9/4Bn38foqKiMGXKFNjZ2aFnz54YP348wsLC8PDhw3Lr+Nh//x184OLigmHDhlU4rWXLluXa2rZtCzU1Ndy5cwfGxsblpq9duxYtWrTApEmTAABeXl4QCoVYv359uW0eOHAgunbtyj1v2rQpgC9/r/r37w8AyMnJwYQJE1BcXIzAwECRQkvoQoJ658mTJ7h16xZmzZqF0NBQ7hEcHIwmTZpwJz6zs7ORl5fHve7OnTto3bo11NTU0LFjR6SmpqJt27bc49GjR9ixY8dnf3g+6NixI+Lj40V+UG7dugV5eXm0adMGenp6SEhIEHlNYmJihctq164dnj59ivz8fK4tNDQUo0aNKjdv+/btUVpaipSUFK4tJiYGlpaW0NbWhry8POLi4rhpjDHEx8ejQ4cOX9ymT/lv7oSEBOjr6wMAOnTogNu3b4tMj4uLq/T6OnbsiCdPnoh8DqdPnxbZU/qcql79165dO9y7d0+kzc7ODhEREV987V9//QVbW1usX78eY8eOhbGxMR4/flztK+jU1NREtvvjR0XFXEZGBg4ODggMDERJSYnItKysLBw+fJjbIwwKCkJkZCTWrl2L5s2bl1uWsrKyyPrU1NQq/b0qLCzExIkTuXukqOCUR0WnngkLC4OCggLc3d2hq6vLPTp37oyRI0fi3LlzyM/PB5/Px5IlS5CcnIxz585h165dcHd3BwC4urri8ePHWLVqFVJTUxEVFYXly5dDUVGxUkXn+++/R1ZWFlauXIlHjx7h4sWLWL9+PRwdHaGiogJnZ2c8fvwY3t7eSEtLQ2BgoMhVWR+zsLBA27ZtsXjxYiQlJSE2Nhbbtm2DpaVluXk7duwICwsL/Pzzz7h79y7u3buHdevWwdzcHMrKynBzc8PatWtx8eJFPHr0CCtXrkRmZibGjBlT7fc7MTERmzdvRlpaGvz9/XH58mXuJtcpU6bgyJEjOHjwINLT03Ho0CEcPXqUe5+/ZNKkSfjrr79w6NAhZGRk4MiRI9i2bRu0tLQq9XolJSW8ffsWKSkpKCkpQWlpKXJzc0UOD35s/Pjx+OuvvxAWFoaMjAxs2rQJL168QM+ePb+4LjU1NSQkJODevXtIS0vD5s2bcenSpU+uSxKmT58OoVAId3d3XLlyBZmZmTh79iw8PDygp6cHZ2dnpKamwtvbG05OTujcuTNyc3NFHh//cfOxynyvXr16hW3btiEzMxPr16+HvLw8t9wPh6UJHV6rd8LDw2FrawtVVdVy01xdXREUFAQ5OTmoqanByMgILi4uUFJSgoeHByZMmAAAaNWqFfbu3YuNGzfCwcEBTZs2ha2tbblLQT+lZcuW2Lt3LzZs2AAHBweoq6vDyckJM2bMAABoa2tj7969WLNmDQICAtCtWze4uLiUOxQDvD98tXPnTqxcuRKjR4+GiooKHB0dMW3atArXvWHDBqxcuRLu7u6Ql5eHjY0Nd5f5vHnzwOPxsGTJEhQVFcHIyAgHDhz4qr9GLSws8PjxYzg4OKBNmzbYvn07d3hn4MCB+OWXX+Dr64s1a9agbdu2WLVqFezt7Su1bGtrayxbtgx+fn5YvXo1tLS08Msvv8DOzq5Srzc3N0fnzp3h6OiITZs2QU1NDe7u7jhw4ECFvQvY2dkhNzcXW7ZswcuXL2FgYABfX1+R8zifMnv2bCxduhRubm5QUFCAiYkJfvrpJ2zfvr3cnoekqKmp4eDBg9i5cyeWL1+OFy9e4JtvvoGNjQ2mTZsGBQUF3L59G+/evcOxY8fKXe4MAGZmZggICCjXXpnvVXh4OE6fPo2SkhKMHTtWZB5DQ0MEBweLb2PrMB6r7v4vqbOCg4Oxfv36St+7QSq2ePFivH37tk7d9b9ixQo4OjqKnLMgpCbR4TVCGojMzEzcv38fnTt3lnYU0oBR0SGkgfj2228REBAg1iveCKkqOrxGCCGkxtCeDiGEkBrToItOWVkZsrKyUFZWJu0ohBDSIDToovPs2TMMGjQIz549k3YUQghpEBp00SGEEFKzqOgQQgipMVR0CCGE1BgqOoQQQmoMFR1CCCE1hooOIYSQGkNFhxBCSI2hokMIIaTGUNEhhBBSY6joANQjASGE1BAqOgBu3bol7QiEENIgUNEBFR1CCKkpVHRARYcQQmoKFR0AT58+RVZWlrRjEEJIvUdF53+io6OlHYEQQuo9Kjr/Q0WHEEIkj4rO/1DRIYQQyaOiA0BZWRkZGRl0XocQQiSMig6Abt26AaC9HUIIkTQqOgBMTU0BUNEhhBBJo6IDoEePHgCAixcvQigUSjkNIYTUX1R0AHTq1Amampp49uwZbt++Le04hBBSb1HRAcDj8TB06FAAQEREhJTTEEJI/UVF538+FJ3Tp09LOQkhhNRfVHT+p1evXlBTU8OjR4+QkpIi7TiEEFIvUdH5Hzk5OVhbWwOgvR1CCJEUKjofofM6hBAiWfWq6Ozfvx9+fn7Vfr2lpSUaN26MuLg4ZGdnizEZIYQQoB4VnV9++QVBQUFftQxFRUUMGDAAAO3tEEKIJNSbotOvXz9Mmzbtq5fz4RDbiRMnvnpZhBBCRNWbojN48GCxLGfIkCFQUFDAtWvX8PTpU7EskxBCyHv1puiIi6qqKgYPHgzGGI4fPy7tOIQQUq9Q0amAk5MTAODYsWNSTkIIIfULFZ0KWFlZQU1NDffv38eDBw+kHYcQQuqNelV0nJyc4OHh8dXLUVBQgJ2dHQAgJCTkq5dHCCHkvVpbdBISEmBubi7S9uDBAzg7O6Nbt24YMWIEEhISJLb+D4fYQkJCaLgDQggRk1pXdBhj+PvvvzFx4kTw+XyuvbS0FNOnT8ewYcNw48YNeHp6wsPDA4WFhRLJYWZmBk1NTTx58gQxMTESWQchhDQ0ta7o+Pj44NChQ+XuuYmJiQGfz8f48eMhLy8POzs76Ojo4NSpUxLJISMjg5EjRwIAgoODJbIOQghpaGpd0XFxcUFwcDC6dOki0p6SkoKOHTuKtHXo0AFJSUkSy/Kh6Jw8eRKlpaUSWw8hhDQUta7otGzZssL2t2/fonHjxiJtioqKKC4ullgWAwMDGBgYID8/HxcuXJDYegghpKGodUXnU5SUlFBSUiLSVlxcDCUlJYmu19HREQAQGhoq0fUQQkhDUGeKTseOHZGWlibSlpqaCh0dHYmu18HBAcD7DkCLiookui5CCKnv6kzR6dWrFxhj2L9/P/h8Pk6ePImHDx9yA69JyrfffouePXvi3bt31PM0IYR8pTpTdBo1agRfX19ERETAzMwMu3btwo4dO6ChoSHxdX+4oIBuFCWEkK/DY4wxaYeQlqysLAwaNAjnzp2Dtrb2J+fLy8uDiYkJGGOIi4tDs2bNajAlIYTUH3VmT0eaNDQ00L9/fwgEAoSHh0s7DiGE1FlUdCrpu+++AwD8/fffUk5CCCF1FxWdShoyZAhUVFRw+/Ztid6QSggh9RkVnUpSVFSEvb09ANrbIYSQ6qKiUwWjR48G8H5wt7KyMimnIYSQuoeKThX06NED7du3R05ODi5fviztOIQQUudQ0akCHo/H7e0cOXJEymkIIaTuoaJTRaNGjQKPx0NERATy8/OlHYcQQuoUKjpVpKWlhb59+6KkpARhYWHSjkMIIXVKtYuOUChESkoKrl+/jtjYWKSnp6OhdG7w8QUFhBBCKk+uqi+4ceMGDhw4gKioKJFel3k8HlRVVdGvXz98//33MDU1FWvQ2mTYsGFYsmQJV2zbtWsn7UiEEFInVHpPJz09He7u7li6dCm0tbWxfft2XLp0CQkJCYiPj8f58+exbt06fPPNN5g3bx7GjRtXbiiC+qJJkyawtbUFQHs7hBBSFZXu8HP06NGYMWMGBgwY8MV5GWM4e/Ysdu/ejaNHj35tRompbIefFbl8+TJcXFzQpk0bREVFgcfjSSglIYTUH5UuOoyxKv+wVuc1Nelrio5AIECvXr2QnZ2NkJAQmJmZSSglIYTUH5U+vPZx8Vi6dCkePnxYbp68vDwMGjSowtfUN7KyslwnoLV5b44QQmqTal29FhwcDBcXF5w4cUKkXSgU4unTp2IJVhd8KDrh4eEoLi6WchpCCKn9qn3J9I8//ohly5Zh7dq1EAqF4sxUZ+jq6qJbt2548+YNdu3aJe04hBBS61W76AwfPhwHDx7Ev//+i/HjxyMvL69eH077FC8vLwCAj48PUlJSpJyGEEJqt2oVnQ/FxcDAgDuf4eTk1CDHmenTpw+cnZ1RWlqKxYsXN5gbZAkhpDqqVXQ+/mHV0NCAv78/Bg8ejKlTp4otWF3y888/Q0NDA9HR0dQRKCGEfEa1is7MmTOhpKTEPZeVlcXPP/+MFStWoEePHmILV1doaGjg119/BQD89ttveP78uXQDEUJILVXp+3Sqcr5CR0en2oFq0tfcp/NfjDG4ubnhwoULGDx4MPbv398gz3ERQsjnVLrvteHDh4PH43E3fH584+fHdYvH4+H+/fviT1rL8Xg8bNiwAYMGDcLZs2fx119/YezYsdKORQghtUqli865c+dEnjPGMGLECOzZsweamppiD1YXaWpqYvXq1Zg1axZ++eUX9O3bF23atJF2LEIIqTUqXXS0tLQqbG/VqtUnpzVEI0eOxOnTp3Hq1CnMmzcPR44cgYwMDVtECCEADeImdjweD+vXr0fz5s0RHR2NgwcPSjsSJyIiAp07d0Z4eLi0oxBCGigqOhKgoaGBlStXAgBWrVqFZ8+eSTkRUFxcDC8vL7x+/RrLli1DYWGhtCMRQhogKjoSMmLECFhbW6OgoAA///yztONg7969yM7OBgDk5uZi+/btUk5ECGmIKn1OJygoqFybQCDA8ePHoa6uLtLu6ur69cnqOB6Ph9WrVyMqKoo7x/Nh4Lea9uLFC67ILFq0CN7e3tizZw/c3Ny++lJxQgipikrfpzNw4MDKLZDHK3elW20lzvt0PmX//v3w8vKCpqYmrly5AgUFBYms53O8vLywf/9+DBw4EAEBAZgxYwZCQ0Ph6OiIHTt21HgeQkjDVemiUx/VRNERCoWwtrbGgwcPsH79eri5uUlkPZ9y+/Zt2Nvbc6O56unpISsrC5aWlnj37h0NQEcIqVGVPqdjY2ODFStW4OzZs3QSugpkZGQwZ84cAMC2bdvA5/NrbN3nz5/H6NGjIRAI4OrqCj09PQCAtrY210/e3LlzUVRUVGOZCCENW6WLzu7du6Gjo4Njx45h4MCBcHFxwfbt23Hr1q0GO55OZdnZ2UFHRwdZWVk4duxYuemxsbHw8PBATEyMWNbHGENgYCDGjx+Pt2/fwsnJCb/99pvIPHPmzEHnzp2Rnp7O9RtHCCGSVq3DawKBAHFxcYiKisLVq1eRnp6OHj16oG/fvujbty/atm0riaxiVxOH1z4IDg7GrFmz0LZtW1y6dAlycu+v4YiMjMSkSZPw7t07qKqqIjw8vNp91zHGcPnyZWzatAmxsbEAgNmzZ2PRokUV9gP38OFDDBs2DCUlJfD394eNjU31N5AQQipBLOd0CgoKEB0djatXryI6OhpnzpwRRzaJq8miU1ZWhgEDBiAtLQ0zZsyAlZUVsrKysHDhQvD5fLRq1QrPnj1Du3btEB4eDnV1ddy9exdpaWmwtbXlihQAvH79Go8fP4a2tjbU1dWRk5ODU6dOITg4GHFxcQAAdXV1LFu2DM7Ozp/N5evri19//RUaGho4duwYdHV1Jfo+EEIaNrqQoIaKDgAcOXIEc+fOLdc+efJkLFy4EN999x3u3LmDzp07g8/nIzk5GcD7Kwf/+OMPKCsr4/z585g5cyZev34NAFBVVUVBQQHX6aqamho8PT0xYcIEKCsrfzGTUCiEq6srLl26BEVFRaxduxajR48Wmae0tBT+/v7Iy8tD7969YWZmhiZNmnzt20EIaYAqVXRSU1ORk5MDIyMjkR+yyMhIWFlZSTSgJNV00REIBPD19cXdu3eRkZGBly9fws3NDZ6enuDxeMjOzsbw4cO5Hgw0NDQgFAqRn58PAwMDWFtbY9u2bWCMQVtbG3l5eXj79i0UFBRgZWUFOzs7DBkypMoF4e3bt1i8eDF3vum7777DzJkzoauri8ePH2PatGmIj4/n5peTk4O2tjaaN2+OFi1aYPLkyejVq5f43ihCSL31xaITFBSEgIAAtG/fHvfv34eXlxesra0BvO/cMiQkpEaCSkJNF53KSE5ORmBgICwsLDBgwABkZWXB3d0dqampAN7fBzV//nzMmTMHPB4PL1++hJKSksigetXBGMOhQ4ewbNkyvHv3DgBgamqKlJQUvHnzBtra2rCzs0NMTAzi4+NFLh4xNTWl/twIIZXDvmD48OGsqKiIMcZYRkYGGzlyJNu3bx9jjDEHB4cvvbxWy8zMZLq6uiwzM1PaUT4rLy+PjR07lnXp0oWdOXNGoutKSkpiCxcuZLq6ukxTU5NpamqyCRMmsFevXnHzFBYWsuTkZHb58mWmra3N2rZty96+fSvRXISQ+uGLezq2trY4deoU97ywsBBz5syBjo4Orl27huPHj0u8MEpKbdzT+RyhUFhjwyS8ffsWJ0+ehKysLEaOHPnJUVCtra2RmJiIY8eOoXfv3jWSjRBSd33xF6x58+YiI4EqKytj9+7dePXqFZKSkiQajoiqyXF5lJSUMHr0aDg5OX122O0PvRmI6x4jQkj99sVfMW9vbzRv3lykTU5ODt7e3ggMDJRYMFI3fCg6N27ckHISQkhd8MWi06pVK7Ro0QIAyl000L17d8mkInVGjx49AAA3b96knikIIV9UpeM1y5Ytg6urKx1WIxwtLS1oaWnh9evX9O+CEPJFVSo6wcHBkJGRgZOTE9atW0cdRRIAQM+ePQHQeR1CyJdVqejo6uoiICAA69evxz///INhw4aJXNlGGqYPRYfO6xBCvqRal0PZ2dnh9OnTcHBwwE8//YQJEyYgLS1N3NlIHUEXExBCKqva1+AqKipi/vz5CA8Ph5ycHOzt7bF582bubnbScOjp6UFFRQWZmZl4+vSptOMQQmqxr77xo127dvD19cXWrVtx4sQJ2NraiiMXqUNkZWW5q9hob4cQ8jlyX56lYpmZmXj48CGSkpK4/2ZnZ+MLHRyQeqpHjx6IjIzE2bNn4eDgIO04hJBaqkpDGyxfvhxJSUlISkpCcXExgPdDHxsYGEBPTw/6+vrQ19eHlpaWxAKLU13rBqc2S0tLw4ABAyAQCBAREQFDQ0NpRyKE1EJV2tN58OAB9PX1YW9vD319fejp6dG4KgQA0L59e4wfPx579+7FihUrcPjw4c92n0MIaZhoEDfa0xGbV69ewcLCAvn5+di/fz8GDRqE4OBgxMTEYNGiReW6UyKENDzVupAgNDQUycnJEAgE4s5D6jB1dXX8+OOPAIBff/0VQ4cOxZw5cxAUFAQvLy8pp5bXg0gAACAASURBVCOE1AbV2tPp27cv8vLyIC8vj06dOsHAwIB76Ovrf/WAYjWF9nTEr7S0FFZWVkhPTwcAtG7dGvn5+SguLkZQUBAGDBgg1XyEEOmq9uG1nJwc3Lt3D4mJiYiPj8e1a9dQVlYGHo+HxMREceeUCCo6knHt2jWsWLECdnZ28PDwwP79+7Fq1Sq0a9cO586dQ+PGjaUdkRAiJWI7p5OTk4N58+bB3t4ezs7O4likxFHRqRl8Ph82NjZISkrCggULMHfuXGlHErugoCA8ePAA48aNg66urrTjEFJriW1UsJYtW2L16tW4cOGCuBZJ6gl5eXmsXbsWALBx40a0adMGWlpasLS0REpKCjdfbm4unJ2d4enpiTdv3kgrbpVdv34dixYtwr59+2BlZYXx48fj3r170o5FSK1UraKTl5dXYXu7du3qzKE1UrN69+6N8ePHAwB3AUpKSgqcnZ2RlpaG3NxcjBkzBleuXEF4eDgcHR2RmZkpxcSV8+7dOyxYsAAAYGpqisaNG+Pff/+Fi4sLSkpKpJyOkNqnWkWnT58+sLS0hKenJ3x8fHD27FmkpqYiNDQUfD5f3BlJPbF69WqkpqYiPT0dSUlJMDc3x7NnzzBmzBiMHj0aSUlJ0NPTQ6dOnfDw4UMMHz4c8fHx0o6N58+fY8+ePSgsLCw3bcuWLUhNTUWnTp1w9OhRXL9+HR06dEBeXh5iY2OlkJaQ2q1a3eAcOXIEiYmJSExMxMWLF7F3716UlpZCVlYWP/30k7gzknpEQUEBwPtDbn/++Se+//577sdZX18fhw8fhry8PKZMmYIrV65gwoQJuHjxIlRUVKSSVyAQYMqUKbhx4wbu3LmDbdu2cdPu3r2LnTt3gsfjYcOGDVBQUICCggKsra2xe/duXLp0CX379pVKbkJqLSYGZWVl7PHjxyw/P18ci6sxmZmZTFdXl2VmZko7SoP1+vVr5ubmxr777juWm5vLtZeWljI7OzumqanJfv75Z6nl8/X1ZZqamtzj33//ZYwx9vLlS2ZpaVlhvsjISKapqcmGDh0qjciE1GqVvnrtzZs3UFVVrVJBe/36NZo2bVqtYlgT6Oq12u3evXsYNmwYGGM4ceIEunbtKpblMsawdetWaGpqYsyYMZ+cLz09HYMGDcK7d+9gY2ODM2fOoHXr1jh58iQmTZqEW7duQU9PD2FhYVBWVuZeV1xcjM6dO4PP5yMhIQEaGhpiyU1IfVDpczpubm7Ys2cPCgoKvjhvXl4edu7cCTc3t68KRxo2Q0NDeHh4QCgUYvHixWLrAePMmTPYsGED5s2b98mhGIRCIRYsWIB3797B0dERe/fuhYmJCbKzszFw4EDcunULWlpaCAoKEik4wPuxpszMzMAYw+XLl8WSmZD6otJ7OkVFRdiyZQuOHj0KU1NT9OvXD506dYK6ujoYY8jLy8ODBw9w48YNxMbGYuTIkZg7d265L2RtQns6tV9RUREGDBiAp0+folevXtDR0UG7du0wbty4ap3nYYzB1tYWCQkJAIAOHTrgzJkzUFRUBADcv38fYWFhOHHiBFJTU9G8eXNERkZCQ0MDDx48wNChQ8Hn86Guro7Q0FDo6OhUuJ4dO3ZgzZo1GDt2LDZu3Fj9N4CQeqbKN4fm5ubi8OHDuHDhAu7fv8/99SknJ4fOnTvD0tISo0aNQsuWLSUSWJyo6NQN//77LyZNmoSysjKuzdLSEkFBQVXuyfr8+fMYN24cWrRoATU1NSQnJ8PT0xMzZszAihUrcPToUW7eZs2aYdu2bbC0tOTaAgICEBgYiHXr1sHExOST67l79y6GDBkCTU1NxMTEUI/bhPzPV/VIIBQKkZ+fDx6PB3V1dXHmqhFUdOqOjIwMPHjwAE+fPoW3tzdev34NPz8/DB06tNLLYIzB3t4et27dwrJly9CrVy/Y29uDMQY1NTW8evUKCgoKGDVqFEaMGAFzc3PIyVVvnEOhUIhu3brh5cuXuHjxInR0dFBYWAglJSXIyIjtnmxC6pxqjxwKADIyMnSSlNSINm3aoE2bNtxzLy8v/PLLL7C0tOQOjX3J5cuXcevWLairq2PcuHFo0qQJpk+fju3bt+PVq1fo06cP1q9fjw4dOnx1XhkZGfTr1w+hoaEICgpCXl4eQkJCYGZmBn9/f6ldAk6ItFWr6Ojr63/ycIG8vDxatmwJOzs7zJo1C7Kysl8VkJD/cnNzQ2BgIO7fv4+dO3di/vz5X3zNmzdv4O3tDQCYMmUKN/jgvHnzICcnh/bt2+O7774T62Gw/v37IzQ0FHv27OHaoqOj4eLigsDAwDp5dICQr1Wtw2tHjhzB1q1bMWvWLHTr1g0AcOfOHfj4+GDMmDHo1KkTtm/fjkGDBtXqzh3p8Frddf36dTg5OUFBQQH9+/fHixcvIBQKYWtrC2dnZ7Ro0YKb98aNG5g1axYyMzPRrFkzXL16tUb2NHJyctCvXz+UlpbC2dkZjo6OmD9/Ph4/fgx9fX0cOnQI33zzTbnXFRQU4NSpU8jOzoanp6dIr9ylpaXg8XiQl5eXeH5CJKI6N/fY2NiwyMjIcu2XLl1iQ4YMYYwxFhsby/r371+dxdcYujm0bps5c6bIjZsfHm3btmUuLi7M3d2djR07lmlra3M3a6akpNRoxqysLPb8+XPu+dOnT1n//v2ZpqYmMzU1ZTExMdy0mJgY5unpyTp06MBty9atW7npJSUlbOjQoUxfX58dO3aswvWVlZWx0NBQ9uDBA5F2gUDA0tPTmVAoFPMWElI11So63bp1q/DLm5SUxIyNjRlj779sXbt2/bp0EkZFp24rKipiR48eZadPn2YxMTEsIiKCjR8/nisyHx5aWlps9erVrKSkRNqRGWOM5ebmMgcHB6apqcnatGnDvL292ZgxY0QyDxs2jGlqajIDAwOup489e/aIzDN9+nSRXkD4fD6bMWMG09TUZN9++y1bvnw5e/36NYuIiGCDBg1impqabOPGjdLabEIYY1XokeBjHh4ekJeXx7p166CmpgYAePXqFZYsWYKSkhL4+/vD398fx48fR2hoqNj3zsSFDq/VT0+fPsXdu3cBADweD+3atUOnTp2knEoUn8/H6tWr4evry7WpqKhgwoQJcHV1hba2NkaPHo2oqCjMnj0bkydPhoWFBV6/fg13d3ccPXoUb9++RfPmzTFt2jS4uLhg4cKFOHXqFBo3bozS0lIIhUIoKCiI9HYtKyuLsLAw7rA4ITWtWkUnOzsbkydPRkZGBrS0tCAUCpGdnY2OHTvCx8cHjx49wuzZs8vd41DbUNEh0hYeHo5du3Zh0KBBmDhxIvdHHADExsbCwcEBSkpKGDp0KIKDg2FhYYG//voLaWlpIj0qyMvLg8/no2nTpggMDESjRo2wdOlS3Lx5E82bN8esWbPw+PFj7Nu3D7q6ujh9+jSN4Eqkotr36QgEAly7dg0PHz6EnJwcOnXqBHNzcwD/P95Obb+cmooOqe3c3d1x7tw5AO8vw46IiEDnzp0BvL/vKDIyEr///jvi4uKgrq6OQ4cOwcjICMD7e4Xi4+Ohp6cHJSUlFBcXw8bGBqmpqZg5cyaWLFkite0iDVe171KTlZWFjIwM5OTkIBQKwefzuTvGNTQ0an3BIaQuWLRoEff/Y8eO5QoO8P7Q4cCBAxEeHo6wsDD8+++/XMEB3hcpExMTKCkpAXjfJ9zmzZshIyODnTt34s6dOzW3IYT8T7Xu08nNzcWMGTOQmJgILS0tMMbw9OlTtG/fHvv370ezZs3EnZOQBqlLly6YNGkSoqKisHDhwgrn4fF46N69e6WW16NHD0yYMAF+fn7w9fWFj4+POOMS8kXVOrz2448/IicnB1u2bOH6WMvJycGCBQvQsmXLOtPBIR1eIw1RRkYGzM3NoaCggJs3b9JNqqRGVevw2uXLl/Hzzz+LdOrZsmVL/PTTT7h06ZLYwhFCxK9NmzawtLRESUmJSAen2dnZtWJ4cFK/VavoKCgoVNhdCI/HE9uYJ4QQyfkw1lVgYCAYY3j27BlsbW0xfPhwOtdDJKpaRadv375Yu3YtXrx4wbW9ePEC69atg4WFhdjCEUIkw9raGt988w1SUlJw5coVTJ06Fc+fP4dQKMSuXbukHY/UY9UqOosWLcLLly9hZWWFIUOGYMiQIbCyskJxcTG8vLzEnZEQImby8vJwdnYGAEyePBmxsbFo1aoV5OTkEB4ejqysLCknJPVVte/TKSsrw+XLl5GSkoLGjRujY8eO6NOnj7jzSRRdSEAasoyMDPTp0weMMSgoKCA4OBh+fn4IDg7G5MmT8euvv0o7IqmHKl10UlJSKr3QTw3hW9tQ0SEN3cSJExEREYFNmzbBxcWFG/G0SZMmuHHjBpo2bSrtiKSeqfR9OsOHDwePx8PHNeq/FxMwxsDj8XD//n3xJSSESIyPjw+ePHkCPT09AO/vC7KwsMCVK1cQGBiIGTNmSDkhqW8qXXQ+dMVR25SVleGnn35CdnY2WrRoAW9vbygoKEg7FiF1grKyMldwPpg2bRquXLkCX19ffPfdd2jVqpWU0pH6qNIXEmhpaVX6UZPOnDmDVq1a4eDBgzA2NkZYWFiNrp+Q+sbS0hI9evRAbm4uXFxcRK5SJeRrVbvvtdri9u3b3GXaffv2RUxMjJQTEVK38Xg8+Pv7Q19fH8nJyXBxccGrV6+kHYvUE3W+6BQWFnLj3Tdp0gRFRUVSTkRI3aehoYG//voLHTt2xP379zF69Gikp6dLOxapB+p80VFWVuYKTVFREVeACCFfp0WLFjh8+DDat2+P+/fvY9iwYfjnn3+kHYvUcVXuZTovLw8XLlxAcnIyioqKoKysDAMDA1haWkJVVVUSGT/L2NgYV69ehbm5Oa5evUojIhIiRq1bt8bJkycxb948/PPPP/Dw8MCCBQswd+5caUcjdVSVbg719/fH1q1bwePxoK2tDRUVFRQWFiItLQ1ycnKYO3cu3N3dJZm3nLKyMixevBhZWVlQV1fH5s2bKz0iIt2nQ0jlMMawZ88erF69GgKBAHv37sWwYcOkHYvUQZXe0zl27Bi2bt2KRYsWwcnJSeSHvaSkBMHBwfD29kbLli0xZMgQiYStiJycXJ0ZSoGQuorH42Hq1Kng8XhYsWIF5s6dCwMDA7Rr107a0UgdU+lzOgEBAViwYAG+//77cnsSCgoKGDt2LGbOnIkDBw58daiEhARu6OsPHjx4AGdnZ3Tr1g0jRoxAQkLCV6+HEFI1kydPhq2tLQoKCjBlyhRER0dj+/btmD59OpYtW4Z9+/YhKioK1exdizQAld7TSU9P/2IP0gMHDoSvr2+1wzDGcPToUaxfv16kvbS0FNOnT4e7uzsCAwNx5swZeHh4IDIyEsrKytVeHyGkang8HjZt2oTExETcu3cPo0aNqnC+H374AWvWrKnhdKQuqPSezrt37774A6+iooLXr19XO4yPjw8OHTqEadOmibTHxMSAz+dj/PjxkJeXh52dHXR0dHDq1Klqr4sQUj2qqqrYs2cPtLS0YGBgAHd3d2zatAnLli2Dq6srFBQU8Oeff4rlqAepf6p09VpFA7eJk4uLC+bMmYPr16+LtKekpKBjx44ibR06dEBSUpJE8xBCKmZoaPjJG7F79uyJH3/8EcuWLYOOjk6d632eSFaVik5wcDCUlJQ+Of1rb8z8ePjrj719+7bceSRFRUUUFxd/1foIIeI3evRoPHjwALt27cKkSZPQq1cvKCsro3379pg1axbk5eWlHZFIUaWLjqamJg4dOvTF+Vq3bv1VgSqipKSEkpISkbbi4uLPFkBCiPQsXboUycnJOHfuHM6cOcO1q6qqYtKkSVJMRqSt0kXn/PnzkszxWR07dsT+/ftF2lJTU+Ho6CidQISQz5KVlcX+/ftx48YNvHr1Cvfv38fGjRuxdetWODs7Q0VFRdoRiZTUiW5wevXqBcYY9u/fDz6fj5MnT+Lhw4ewtraWdjRCyCfIyMigV69eGDp0KH788UeYmZkhLy8Pu3btknY0IkVVKjqFhYUIDQ1FYWEhgPeXOPv6+sLd3R0//vijxO6dadSoEXx9fREREQEzMzPs2rULO3bsgIaGhkTWRwgRLx6Ph6VLlwIAdu/ejefPn6O4uBhHjhxBQEAA9WLdgFS6G5yMjAy4urqiqKgIx48fx7fffos1a9YgICAAw4YNg7KyMsLDw+Hn5wdTU1NJ5xYL6gaHkJr1YXhsU1NTPH78GC9fvgQANG7cGHZ2dhg3bhx69Ogh8StlifRUuujMnz8fRUVF2LJlCxo3boxXr16hX79+sLa2xubNmwEAfn5+uHr1Kvbt2yfR0OJCRYeQmpWUlIRBgwZBKBQCeN9hr5qaGi5dusTNY2JiwvV8QFe61T+VPrwWHR2NGTNmcJcuX7lyBQKBQORkft++fREfHy/+lISQekFXVxcrVqyAvb09jhw5glOnTuHQoUOIiorCjBkzoKamhri4OEyfPh3Tp0+XdlwiAZUuOgUFBWjWrBn3/Nq1a5CTk4OZmRnX1qRJE+4vGEIIqcjEiRPxxx9/oG/fvtxhtLZt22Lp0qWIjY3FunXr0KRJE5w6dQq3bt2ScloibpUuOlpaWkhJSQEACAQCXLx4ET179oSioiI3z7Vr1+gwFSGk2hQVFTFu3DhMmDABALBlyxYpJyLiVumi4+TkhFWrViEsLAzLli3DixcvMHbsWG56bGwstmzZAjs7O4kEJYQ0HFOmTIGioiLOnTuHO3fuSDsOEaNKFx0PDw8MHjwYq1evRmRkJObNm8fdJ7Ny5Uq4ubmhe/fu8PDwkFhYQkjD0KxZM4wbNw4AsHXrVimnIeJUpZFDP+Xhw4cQCoUwMDAQR6YaQ1evEVJ75eTkwNzcHCUlJThw4ADat28PFRUVNG/enC6prsMq3Q1OcXEx1qxZg4iICDRq1AiDBw/GggULoKysDD09PUlmJIQ0QC1btoSrqyv27dsHd3d3rv3bb7+FlZUVrK2tYWVlRQWojql00fHx8UFkZCQmTZoEGRkZBAUFIT8/n070EUIkZvbs2UhLS8OzZ89QXFyMvLw8ZGZm4sCBAzhw4ACcnZ2xceNGyMjUiR69CKpQdCIiIrBx40b07t0bAGBmZgZXV1fw+Xy6gYsQIhEtWrRAYGAg91wgECAhIQHnz5/HH3/8gcOHD0NBQQFr1qyhPZ46otJ/HuTk5KBDhw7ccyMjIzDGuG4sCCFE0mRlZWFiYoL58+fD398fCgoKOHDgAH799VeI4fQ0qQGVLjoCgQCysrLccx6PB3l5efD5fIkEI4SQz+nXrx/27t0LeXl57N27Fx4eHvRHcB1QpZFDCSGkNhk4cCD27t2LWbNmISIiAnFxcZgzZw4yMjJw8+ZNFBUVwcjICF27doW1tTW0tLSkHbnBq/Ql0/r6+pg/f77IaJ3r16/HlClToK6uLjKvq6ureFNKCF0yTUj98OTJE8yePRvXrl375DwfesLX1dWtwWTkvypddAYOHFi5BfJ4OHfu3FeFqilUdAipPwQCAfz8/HD9+nUYGBige/fuUFZWxp07d3Dy5Elcu3YNHTp0wMmTJ6GqqirtuA2WWG4Orauo6BDSMLx9+xb29va4f/8+bGxs4OfnR5dZSwm964SQek9JSQl79+5F06ZNcebMGaxevRoCgUDasRokKjqEkAahXbt22LFjB3g8Hnbt2gV7e3s8ePBA2rEaHCo6hJAGw8rKCgEBAWjdujVu376NoUOHYvny5UhPT5d2tAaDig4hpEGxsrJCZGQk3N3dwefz4efnBwsLC/zwww+Ii4uTdrx6j4oOIaTBUVFRwdq1axEREYExY8agUaNGOHv2LEaMGIFFixYhLy9P2hHrLSo6hJAGq0uXLti8eTNu3LiB6dOnQ1ZWFkFBQejXrx8uX74s7Xj1EhUdQkiD16xZM3h5eeHs2bOwsLBAfn4+pk6dirS0NGlHq3eo6BBCyP906tQJhw4dwtChQ/H69Wt4eHigsLBQ2rHqFSo6hBDyERkZGWzZsgU6Ojp4+PAh5s6di7KyMmnHqjeo6BBCyH+oqKhg3759UFFRwalTp9ClSxdMmjQJQUFBKC4ulna8Oo2KDiGEVKBjx47Yu3cvdHR0UFBQgNOnT2PRokWwsLBAUFAQ7f1UE/W9Rn2vEUK+IDMzExcvXkRAQADu3r0L4H0PB+PGjcOYMWOgoaEh5YR1BxUdKjqEkEoSCoUIDw+Ht7c314tBo0aN4OjoiF9//RVNmzaVbsA6gA6vEUJIJcnIyMDBwQEXL16Ev78/Bg0ahLKyMhw5cgRDhw5FQkKCtCPWelR0CCGkiuTk5GBjY4MDBw7g4sWLMDIyQkZGBhwcHLBx40bcvn2berH+BCo6hBDyFTp06IDQ0FD88MMPKC0txebNm2FnZwdDQ0N4enrizJkzKC0tlXbMWoPO6dA5HUKImFy4cAEnTpxAVFQUHj9+zLWrqanBxcUFU6ZMQcuWLaWYUPqo6FDRIYRIQEZGBsLCwhAcHIyHDx8CABQUFODs7IyZM2dCS0tLygmlg4oOFR1CiAQxxhAfH48dO3bg1KlTAN4Xnx9++AGzZs1qcJdb0zkdQgiRIB6Ph27dusHX1xeRkZGwt7dHSUkJ9uzZgz59+iA4OFjaEWsUFR1CCKkhurq6+OOPPxAREYEBAwagoKAAs2bNwqZNm9BQDjpR0SGEkBrWpUsXBAUFYeXKlZCRkcHvv/+OmTNnIiUlpd4XHyo6hBAiJRMnToS/vz+aNGmC0NBQWFpaonfv3li+fHm9HVKBig4hhEjR4MGDERYWBkdHR2hoaCArKwt+fn6YPHlyvby/h4oOIYRImb6+Pnbs2IH4+HiEhoaiefPmuHTpEubOnQuhUCjteGJFRYcQQmoJGRkZ9OzZE4GBgdwhtyVLliAjI6PenOuhokMIIbWMkZER9u7dC3l5eQQGBsLc3By9evXC0qVL8eTJE2nH+ypUdAghpBbq378/AgMDMXToUKipqeHJkyf4888/YWFhAS8vL+Tk5Eg7YrVQ0SGEkFrKwsICfn5+uHPnDk6fPo2RI0eCz+dj//79sLW1RW5ubrnXlJWVISwsDNOnT0dsbKwUUn8eFR1CCKnlZGRkYGxsjO3bt+Ps2bMwMTHBs2fPMGfOHO5Cg9LSUuzatQvm5uaYNm0ajh8/Dk9Pz1p36TUVHUIIqUP09fXh6+sLdXV1XLx4EX/88QcyMzPh5OSElStX4unTp+jYsSN0dHSQnZ0Nb29vaUcWQUWHEELqmNatW2PLli0AgPXr12PIkCGIi4uDpqYm9u/fjwsXLmDHjh2QlZWFv78/4uPjpZz4/1HRIYSQOmjw4MGYMmUKBAIBXr9+DWtra5w5cwbW1taQkZFBly5dMHnyZAiFQixcuBAvX74En8//7DILCgqwZs0aifYFJyeRpRJCCJG4JUuWoFGjRtDW1oabmxt4PJ7I9Pnz5+PEiRO4d+8ejI2NAbwfUG7cuHGYOnUq1NXVuXkvXLiAhQsX4unTpwCAAQMGoHv37mLPTOPp0Hg6hJB67MaNG1iyZAmys7Px5s0b7sIDZWVl2Nvbo7CwEJmZmYiLi+PaCwsL4eLigk2bNok9DxUdKjqEkAaCMYabN2/i999/x8WLF0WmKSgoYMGCBbC2tsaAAQOgqKiIW7duQVVVVawZ6PAaIYQ0EDweDz169MDBgwdx48YNxMbG4ptvvoG2tjY6derEjWJqbm6O6OhohISE4IcffhBrBio6hBDSAPXs2RM9e/ascJqbmxuio6MRFBQEd3f3cueKvgZdvUYIIUTEh6537t27h4SEBLEum4oOIYQQEY0bN8bo0aMBAAEBASLTBAIBoqKiqr1sKjqEEELKcXV1BQAcOnQIM2fORF5eHu7cuQMHBweuIFUHFR1CCCHldOrUCevWrUPjxo0REhKCfv36wdbWFnFxcWjVqlW1l0tFhxBCSIXGjRuHs2fPwtzcHPn5+QCAyZMnl7vcuiro6jVCCCGf1L59exw5cgRnz55F27Ztoaen91XLo6JDCCHks2RkZGBjYyOeZYllKYQQQkglUNEhhBBSY6joEEIIqTFUdAghhNQYKjqEEEJqDBUdQgghNaZBXzItEAgAAM+ePZNyEkIIqXtatWoFObmqlZEGXXRyc3MB/H8fQ4QQQiqvOgNgNuiRQ9+9e4e7d++iRYsWkJWVlXYcQgipU6qzp9Ogiw4hhJCaRRcSEEIIqTFUdAghhNQYKjqEEEJqDBUdQgghNYaKDiGEkBpDRYcQQkiNoaJDCCGkxlDRqWOOHj0KQ0NDmJiYcI+QkBBpxxKbhIQEmJubc89LS0uxbNkymJmZoXfv3ti9e7cU04lPRdvZpUsXkc914sSJUkxYfVevXoWTkxNMTU1hbW2Nv/76CwDw5s0bzJ49G927d0e/fv1w7NgxKSf9Op/azszMTOjr64t8ll5eXlJOW32RkZEYMWIETExMMHjwYG47q/3dZKROWbFiBduwYYO0Y4idUChkR44cYd27d2fdu3fn2jdu3MhcXV1Zfn4+y8zMZEOGDGEhISFSTPp1PrWdCQkJrE+fPlJMJh5Pnz5lJiYm7MyZM0wgELD4+HjWs2dPdunSJTZ37lw2Z84c9vbtW5aYmMjMzc3Z9evXpR25Wj63nadPn2ajR4+WdkSxyMnJYV26dGEXLlxgjDF29+5dZmRkxO7evVvt7ybt6dQx9+7dg4GBgbRjiJ2Pjw8OHTqEadOmibSHhITA09MTTZs2hba2Njw8PLi/tOqiT23nvXv3oK+vL6VU4vPkyRMMHz4c1tbWkJGRgbGxMczMzBAdHY2IiAjMmTMHioqKMDAwwOjRo3H48GFpR66WT23nrVu36s1nCQDffPMNoqOjYWlpCaFQiPz8fMjKyqJJkybV/m5S0alDBAIBHj58iOPHw6zY3AAACiZJREFUj8PCwgLW1tbYs2cPWD3oycjFxQXBwcHo0qUL1/bmzRvk5uZCR0eHa2vfvj2SkpKkEVEsKtpOAEhMTEReXh5GjBiBPn36YPbs2cjJyZFSyurr0aMHfvvtN+55fn4+YmNj8c0334DH46Fdu3bctA4dOtTZz/JT29m5c2ckJiYiOTkZQ4YMgYWFBby8vPDmzRsppv06ysrKKC4uhpGRESZOnAhXV1doaGhU+7tJRacOycvLQ5cuXeDo6Ijz589zfzUfPHhQ2tG+WsuWLcu1vX37FgDQuHFjrk1RURHv3r2rsVziVtF2Au+3y9TUFH/++Sf++ecfNG7cGDNmzKjhdOJVUFCAadOmoWvXrjA0NESjRo3A4/G46XX9s/zg4+0cNGgQVFVV0bdvXxw9ehQhISHIzs7GsmXLpB3zqygoKCAuLg5Hjx7FsWPHcODAAQDV+2426KEN6poWLVogMDCQe25gYAA3NzecOXOmXg7PoKioCAAoKSnh2oqLi6GkpCStSBKzZMkSkeeLFy+Gubk5srOz0bp1aymlqr60tDRMnz4dOjo62LhxIx49eoTS0lIwxrjCUx8+y/9up4yMDDZv3sxNV1FRwdy5c/H999+jrKysyj0y1xYyMjJo1KgRjIyMMGbMGNy9exdA9b6btKdThyQnJ8PHx0ekjc/nQ0FBQUqJJKtp06Zo0aIFUlNTuba0tDSRXfr6YuvWrXj06BH3nM/nA0Cd/Gxv3LiBMWPGYPDgwfDx8YGCggLatm0LxhiysrK4+VJTU+v0Z1nRdhYXF8Pb2xsvXrzg5uPz+ZCTk6uTw6fExMTAyclJpK20tBSqqqrV/m5S0alDVFVV4e/vjyNHjkAoFOLu3bsICAgo94+iPrG3t8eOHTuQl5eHrKws+Pn5wd7eXtqxxO7hw4dYt24d3rx5gzdv3mD16tUYMGAANDQ0pB2tSjIyMjB16lTMnj0b8+fP5/ZqmjRpgsGDB2PTpk0oKirCgwcP8Pfff9fZz/JT26moqIioqCj8/vvvePfuHZ4/f45NmzZh5MiRIocW6woDAwPk5OTA398fAoEAt27dwrFjxzBq1KhqfzdpPJ06Jjo6Ghs2bEBaWhrU1dXh4eFRrw6tXb9+HTNmzEBsbCyA97vv69at+7/27i+kqT6MA/gXsyGcjZxj1bCxVReiJbRWTELctBr9gYggymAz2nJZMSQRLVJnjXI1pBZbGgWiFSE1kawoR1EuZEFlFwkldNMfsGIzYqFrc++FdNjUnL3tna/z+cC52M45z+/ZgbNn53fOfj/cv38fo6Oj2LVrF8rLy2flCRxp/Of0+Xwwm81wu90IhUJQKpWora3FggULZjjTP3P69Gm0tLRM6GbZs2cPSktLUV9fj6dPn4LD4WD//v3QarUzlOnfmepzFhcXw2w24/nz50hJScHWrVtRVVU1K69agbEnK81mM96+fQuRSASj0Qi1Wv2vz00qOoQQQhKGutcIIYQkDBUdQgghCUNFhxBCSMJQ0SGEEJIwVHQIIYQkDBUdQgghCUNFhyS1oqIiZGVlTboEg8GZTm9GeTwe9licOHECHz58QFZW1qSDNjqdTigUipgxrVYrG/PRo0f/RdpklpudAwER8geOHDky6agNs3UcrHhzuVzg8/kYGhr661hlZWUoKSlBfn5+HDIjyYjOOpL0GIaBUCic6TT+tzIyMsAwTFyKDsMwYBgmDlmRZEXda2ROu3DhAvR6PXQ6HeRyOTv195UrV1BYWAiZTIbi4mL09fWx+wSDQVgsFuTl5UGhUODSpUvYsWMHnE4nAECj0cBisUS1M767aar41dXVMJlMOHr0KGQyGYqKiuBwONj1oVAIdrsdKpUKMpkMGo0GAwMD+Pz5M3JyctDb2xvV9saNG+M6NbTT6Zy0u1Kj0cStDZK8qOiQOa+npwdr1qxBe3s7lEolbty4gdbWVtTV1aGjowNKpRIlJSV4//49gLFC1dnZCYvFgpaWFjx8+BADAwPTbi9WfAC4efMmFi1axA6ueP78eXY4ebvdjra2Nhw7dgwdHR0QiUQwGAwQCATIy8tDV1cXG6evrw+Dg4NQq9VxOlrAli1b4Ha72eXcuXNITU2FXq+PWxskeVHRIUmvoaEBMpksavF4POz6tLQ0GAwGLF++HBkZGWhubkZFRQVUKhWkUikOHDgAuVyO69evIxwOo729HWVlZVAqlcjOzobVakUoFJp2PlPF/0UsFqO8vBzLli3DwYMHkZ6ejtevXyMcDrPTXavVakilUtTW1kKtVuPbt2/Ytm0buru72akRurq6oFKpwOPxpp3fzp07Jxyvurq6qOMlFAohFAoRDAZx8uRJHDp0CEqlctptkLmL7umQpGcwGCYMuR45g+eSJUuQkjL2+8vv9+PTp0+oqamJ+qINBALgcDjwer3wer1YsWJF1P6LFy+eVi6x4v8ikUii9mMYBsFgED6fD16vF7m5uew6LpeL6upqAGNdaSaTCW63GwUFBbh37x7q6+unldsvNpstalppAHjw4AEuX74c9d7IyAgOHz6MVatWoays7I/aIHMXFR2S9Ph8/oQv8UiRQ86Pjo4CGLs6ysnJidouLS0tanreSJEFY7zIR7NjxZ8qXjgcxvz58wHgt8PHMwyD9evX4+7du+BwOAgEAigoKPhtbpMRiUQTjpdAIJiwXU1NDfx+P86cOTPrp5ogiUPda4RE4PF4EAqFGBwchEQiYZerV6+ip6cHDMMgMzMz6sa/z+fDx48f2dccDgffv39nX0feq4kVfzr5CQQC9Pf3s+8NDw9j3bp1ePXqFYCxie+ePHkCl8uFTZs2TVkQ/62Wlha4XC7Y7XZwudy4xyfJi650CBlHr9fD4XBg4cKFyM3Nxe3bt3Ht2jW0tbUBAEpLS9HY2AixWIylS5fCarUiEAiw+69cuRKtra3Yvn07uFwuGhoaor74Y8WPZe/evXA4HMjMzIRUKkVTUxN4PB6ys7MBAPn5+UhNTcWtW7cmdInFQ29vL6xWK06dOoX09HR8+fKFXUePppNYqOgQMo5Wq8Xw8DDOnj2Lr1+/QiqVwmazQS6XAwB2796NoaEhmEwmjIyMQKfT4eXLl+z++/btw7t376DT6cDn82E0GqOuhGLFj0Wn08Hv9+P48ePw+/1YvXo1mpub2cI2b948bN68Gd3d3Vi7dm0cj8yYzs5O/Pz5E5WVlRPWvXnzJu7tkeRCM4cSEgcKhQJVVVWTjnwwE4xGIyQSCSoqKn67jcfjgVarxYsXL+L+h86srCw0NTWhsLAwrnHJ7EdXOoQkkWfPnqG/vx+PHz+O+r/OVLxeL8LhcFzuzfj9fvz48eOv45DkRQ8SEJJE7ty5A5vNhsrKSojF4mnts2HDBjQ2Nsal/YsXL9K4a2RK1L1GCCEkYehKhxBCSMJQ0SGEEJIwVHQIIYQkDBUdQgghCUNFhxBCSMJQ0SGEEJIw/wBW697No7aoPQAAAABJRU5ErkJggg==\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+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd1zV1f8H8NflMmWryBIVB4IICiIOXAy1FHFk5M7QcqXmyFzUN9NwpCWOHJGmZuYANU1xpJa4QAsElI2CIiCXzWXde35/+LufvAF6uVy49+L7+XjcR93PfH8A7+uecz6DxxhjIIQQQhREQ9kFEEIIaV4oWAghhCgUBQshhBCFomAhhBCiUBQshBBCFIqChRBCiEJRsKgRoVCI7du346233oKTkxMGDhyI5cuXIz09XSHb37ZtG8aNGwcACA0NRZ8+fQAAmZmZ6Nq1KxITE2XajkAgwG+//aaQmkj9MMZw9OhRVFRUNHhbt2/fRteuXVFaWqqAyhpPdXU1QkJCMGrUKPTo0QODBw/GqlWrkJ2dzS2zfPlydO3atdbX1KlTX7uPJ0+ewN7eHjNmzKh1/r179zBx4kS4uLjAx8cHe/bsgVgsVtgxqhsKFjUhFAoxbdo0nDt3DosXL8b58+exbds2VFZWYvz48fj7778bvI+AgACEhIQ0eDubNm3CxYsXG7wdUn+RkZEIDAxEdXW1sktpEtXV1QgICMDx48fx8ccf48yZM9i8eTPS0tIwadIk5OXlAQBWrVqF69evS722b98ODQ0NeHp6vnY/p0+fRvv27XHjxg1kZWVJzXv69ClmzpyJXr16ISwsDCtWrEBISAgOHDjQKMesDihY1MS2bdtQUFCAX3/9FcOGDYO1tTV69uyJLVu2YNiwYVi2bBmqqqoatA99fX2Ympo2uFa65lZ53rSf/b59+5CUlISff/4Zw4cPh42NDdzc3LBnzx4IhULui5KhoSHMzMy4F5/Px5o1a9CvXz988MEHr93Pb7/9Bn9/f5ibmyM0NFRq3rlz52Bubo6lS5eiQ4cO8Pb2xvTp03Hy5MlGOWZ1QMGiBsRiMY4dO4Zp06bBwMCgxvyFCxciIyMD169fBwBERUVh/PjxcHZ2xoABA7Bx40aIRCIAgEgkwo4dOzBkyBC4uLhg6tSpSEpKAiDdFfYqz58/x5IlS9C3b190794dPj4+OHbsGLeNsLAwhIeHo2vXrgCAyspKBAcHw8vLC05OTpg4cSKio6O57Xl5eWHjxo0YMmQIBg8ejC+++AKTJk2S2ufJkyfh5eVV6wdneXk5vv76a3h4eKBXr16YM2cO1w3CGMOBAwcwfPhwODk5YfTo0bh27Rq37vLly7F27VqsWLECPXv2xKBBg/D777/jwoUL8PHxgYuLCxYvXozKykru+GbPno21a9fCxcUFgwYNwv79+6XquX79Ot5991307NkTnp6e+OGHH7i6Q0NDMW7cOOzduxcDBgxAnz59sHTpUpSVlXHrX7t2DaNHj4azszNGjhyJEydOcPNetX5mZiamTZsGAHB1da3xAViX8+fPw8/Pj9vfpUuXpOafPn0aXl5ecHZ2xgcffCD1jf2vv/6Cv78/nJ2d0aNHD0ybNg0pKSkA/u1CPX/+PN5++23u7y01NbXWOkJDQ+vsrqrrWI4fP4533nkHLVu2lJpuYGCAPXv24P333691vZUrV6K6uhrr168Hj8d75c8nNjYWKSkp6Nu3L3x8fBAWFib1dzh06FAEBQVJrcPj8VBeXv7K7TZrjKi85ORkZmdnx6Kjo+tcZvjw4WzLli2surqaubu7s23btrGMjAwWERHB3Nzc2NGjRxljjG3dupX16dOHhYeHs7S0NPbpp58yT09PVl1dzYKDg9nYsWMZY4ydOHGCubu7M8YYy8jIYHZ2diwhIYExxlhAQACbPn06e/DgAUtPT2fr1q1j3bp1Yzk5OaykpIQtXLiQzZo1i+Xk5DDGGFu1ahUbNGgQu3r1KktOTmaBgYHMxcWFZWdnM8YY8/T0ZG5ubiw6OprFxMSwu3fvsq5du7KnT59yxxcQEMA2b95c67EvW7aMeXl5sYiICJacnMwCAgLYhAkTGGOM7dy5k7m5ubEzZ86w1NRUFhwczBwcHNiDBw8YY4x99tlnzNHRke3cuZM9fvyYrVixgvXs2ZP5+/uz2NhYduXKFebs7MyOHTvGGGMsODiYOTo6srlz57LExEQWFhbGnJycWGhoKGOMsTt37jAHBwe2Y8cOlpqayk6dOsVcXFzYoUOHuJ+ro6MjmzVrFktMTGQXLlxgzs7OLCQkhDHGWGJiInN2dmaHDx9mjx49YmfPnmW9e/dmZ86cee361dXVLDw8nNnZ2bHHjx8zoVD42r+tGzduMHt7e/bDDz+w9PR0tn//fubo6MiSkpLYrVu3mJ2dHfP392f3799nMTExzMfHhy1atIgxxlhmZiZzdHRkISEh7PHjx+zvv/9mfn5+bNasWVJ/NyNHjmR37txh0dHRbPjw4dz8/xIKhSwnJ6fWV23HIhQKmZ2dHTt//vxrj/NlBw8eZHZ2duzKlSsyLb9u3To2ZMgQxhjjfiY3b96sc/nS0lLm6enJVq1aVa+6mhMKFjVw9+5dZmdnxx49elTnMu+99x4LDAxk+fn5zM7Ojv3666/cvOjoaJaZmcnEYjHr27cv279/PzevuLiYBQUFsby8PJmD5aeffmKPHz/mtiEQCJidnR2LjIxkjL34sJ4/fz5jjLHCwkLm4ODAzp07xy0vEonYiBEj2JYtWxhjL4IlMDBQ6ni8vb3ZDz/8wBhjLC8vj3Xr1o0lJibWOO6ioiLWrVs3dvHiRW7a48eP2aZNm1hFRQXr06cPtx2JGTNmsMWLF3O1jhgxQupnZWdnxyIiIrhpU6ZMYUFBQYyxF8HSq1cvVlxczM1fu3YtGz9+PGOMsfnz57OPPvpIan+7d+9mAwYM4H6udnZ2XOgyxti8efO4epYtW1bjA+n777+X+r28an3JB19JSUmNn1Vt5s+fz+bOnSs1befOnSwmJobbVmxsLDdvx44d7O2332aMMfbo0SN24MABqXVDQkKYt7c3Y+zfvxtJKDL24m/Hw8NDptpe59mzZzV+V6+TmJjInJyc2FdffSXT8tXV1czDw4OtXbuWe9+3b1+2dOnSWpevrKxks2fPZu7u7iwrK0vmupob6gpTA8bGxgCAkpKSOpcpLi6GqakpTExMMGvWLAQGBmLgwIFYvXo1SktLYW1tjfz8fAgEAjg5OXHrGRgYYPny5TW6El5l4sSJ+Oeff7BmzRoEBATg7bffBgCuu+1laWlpEIlEcHFx4aZpaGjAxcWF64IDgHbt2kmt5+fnh7NnzwIAfv/9d3Tu3BldunSpdfvV1dVSx2RjY4OlS5eiqKgI+fn56Nmzp9Q6vXr1QnJycq371tXV5bYhoa2tzXWFAYCDg4NUl6SzszN3LElJSbXuLycnB0VFRQBejGWZmZlx8w0MDLjxsaSkJJw8eRIuLi7ca8eOHUhLS+OWf9X69ZWSkiL1swOAOXPmSE17+edjZGTEdfG0a9cOw4YNw549e7Bs2TKMHz8eW7ZsqXE2VIcOHaRqrevEgtOnT0sd98uv06dP11heMh5YWFgo07FWVFRg8eLF6NChAz799FOpeVFRUVL7+/zzzwEAN27cQG5uLoYNGwYA4PP58Pb2xoULF1BcXCy1jfLycsyfPx+3bt3Crl27YGFhIVNdzZGmsgsgr9e+fXu0bNkSd+/eRbdu3WrMFwgESEtLQ48ePQAAixcvxtixY3H58mVcu3YNAQEBWLhwISZPngwAr+1TfhXGGGbMmIGcnByMHDkS/fr1Q+fOnfHWW2/VuryOjk6d23n5A0jygS7h5+eHHTt24PHjxzhz5gxGjRpV63a0tbUB1H5M/91mXfvW0tKqscyrfkZ8Pl/qvVgshoaGRp37ZP/fHy/ZZ237kxCJRJg6dSomTJhQ5zKvWr++tLS0Xvv3IDm2/0pMTMR7770HDw8P9O7dG++++y6io6Nx+PDhV9bL6jjBwMvLi/sb/q9WrVrVmKatrQ0HBwfcv3+f+3Lzsr179yIvLw/Lly8H8OJsxcePH+PEiRM1/i67d+8uNdgu+eIgCbSXx2okfz9nzpzBxIkTAbz40jdr1iwkJiZi3759Nb5cvGmoxaIGNDU1MWHCBPz4448oKCioMX/r1q2wtLTEwIED8fTpU3zxxRewtrbGzJkzcfDgQXz44Yc4ffo0DA0N0apVK8THx3PrlpeXo3///lKD6a+SnJyM27dvY+/evZg/fz6GDh3KfXOTfGC8/EHVvn17aGlpSZ0OzRhDdHQ0OnbsWOd+OnToAGdnZ5w4cQL379+Hr69vrcu1bdsWfD5f6piePHkCd3d3lJaWok2bNvjnn3+k1vn7779fue/XSUpKkmrB3L9/H/b29gCAjh071rq/Vq1acS3PV+nUqRMePXqE9u3bc687d+7U+LCuS32/NHTo0AFxcXFS02bMmFHjhITahIWFoVu3bti+fTvef/999O7dG5mZmXKfmWZgYCB13C+/ajtpBQDGjh2LEydOQCAQSE0vKirCgQMHuFquXbuGgwcPYvny5ejcuXON7ejq6krtr1WrVhAKhbh06RL8/f1x8uRJ7nXq1Cm0bduWO6lCJBJh7ty5SElJwcGDB9/4UAEoWNTGnDlzYG1tjYkTJ+LixYt4+vQpYmJisGzZMpw5cwabNm2ClpYWWrZsifDwcKxbtw6PHj1CXFwcIiIiuK6N6dOnY+fOnbh69SrS09Pxv//9D4aGhnBwcJCpDiMjI/D5fJw9exZPnjxBREQEPvvsMwDgPmxbtGiBJ0+e4MmTJ9DT08OUKVMQFBSEa9euISUlBV999RUyMjLg7+//yn2NHj0a+/btg6ura53dCgYGBhg/fjyCgoIQGRmJ5ORkfPHFF+jatSvMzc3x0UcfYdeuXTh79izS09Oxc+dOXL9+XaaL4ury/PlzfPnll0hNTUVYWBiOHj3KfaOdOXMm/vrrL+zcuRPp6en4/fffsWfPHkydOlWmD/2AgABcvXoVu3btwqNHjxAeHo6vv/661m/stWnRogUAIC4uDqWlpRCJRMjNza3zDKX3338fly5dwqFDh/D48WMcOHAAkZGRGDhw4Gv3ZWpqitTUVERFRSEjIwP79+/HsWPHpEK3sU2aNAm2traYPHkyLl68iIyMDERERCAgIAC6urqYPXs2BAIBVqxYgQEDBsDHxwe5ublSr/+GksTFixdRVlaGDz74AHZ2dlKvyZMn4/79+0hMTMTPP/+MO3fuYO3atWjVqhW3Xck1NG8i6gpTE9ra2vjxxx+xf/9+fPfdd8jMzISRkRE8PDxw4sQJrh9bV1cXu3fvxvr16zFmzBhoaWnB29sbK1euBPDi22hpaSk39uLq6ordu3dzXUqvY25ujjVr1mDnzp3YuXMnrK2tMWnSJBw7dgxxcXEYNGgQxo4di0uXLmHEiBG4dOkSFi9eDB6PhxUrVqC0tBROTk44cOCAVN97bUaMGIF169bBz8/vlcstX74c69evx8cffwyRSAQPDw+sX78eADBlyhSUlZVh06ZNyMvLg52dHXbt2gU3NzeZjrc2Xbt2hba2NsaNG4fWrVvjiy++wPDhwwG8GH/Ztm0btm7dip07d8LCwgJz587F9OnTZdp29+7dERwcjODgYGzfvh1mZmb46KOP8OGHH8q0vp2dHTw9PREQEIAlS5Zg6NCh8Pb2RlBQUK2nkru4uGDDhg3YsWMHNmzYgI4dO2LHjh3o1KkTnj9//sp9TZ06FQkJCZg9ezZ4PB66deuG//3vfwgMDMSzZ89kqrehtLS0EBISgj179uCbb77Bs2fPYGpqioEDB2L+/PkwNTXF7du3kZeXh+vXr2PAgAE1tmFtbY0//vijxvTffvsN7u7utbZu33nnHQQHB+PEiROIiYkBYwzz5s2TWsbExAS3b99W3MGqER6Tt91KSCNLSkrCO++8g+vXr8PIyEjZ5QB4cR3LlStXZL5GRBX88MMPsLKywogRI5RdCnlDUIuFqJy8vDxERkbi0KFD8PX1VZlQUUclJSX4/fffFXKrHkJkRWMsROWUlZVh5cqVKC8vx+LFi5VdjlozMDDAkSNHFHKrHkJkRV1hhBBCFIpaLIQQQhSq2QdLdXU1MjMz35jbiBNCiLI1+2B59uwZvL29m+z0R0IIedM1+2AhhBDStChYCCGEKBQFCyGEEIWiYCGEEKJQFCyEEEIUioKFEEKIQlGwEEIIUSgKFqKS0tLS6nyGCCFEtVGwEJXz8OFDDBgwgHuGDCFEvVCwEJWTlJQEAEhJSVFyJYQQeVCwEJVTUFAA4MVzywkh6oeChagcChZC1BsFC1E5kmApLi5WciWEEHlQsBCVIwmW0tJSetwBIWqIgoWoHEmwANRqIUQdUbAQlUPBQoh6o2AhKuflYKEBfELUDwULUTkULISoNwoWonKoK4wQ9UbBQlSKUCiUukcYtVgIUT8ULESlvNxaAShYCFFHFCxEpRQWFkq9p2AhRP1QsBCV8t8WC42xEKJ+KFiISqFgIUT9qUywREREYNy4cXB1dcXQoUNx5MgRAEBGRgbs7e3h4uLCvVatWqXkakljkQSLrq4ugJpdY4QQ1aep7AIAICsrC/Pnz8eGDRvg7e2N2NhYzJw5E9bW1igtLYWzszOOHj2q7DJJE5AES7t27ZCYmEgtFkLUkEoEy5MnT+Dr64uhQ4cCAJydneHu7o579+6huroa9vb2Sq6QNJX8/HwAgI2NDQULIWpKJYLFzc0Nbm5u3PuCggJERUVh9OjROHLkCMrKyjB8+HCUlpZi8ODB+Oyzz2BkZKTEikljkbRY2rdvD4C6wghRRyozxiJRXFyMOXPmoEePHvD29oaRkRE8PDxw/PhxhIWFISsrC4GBgcoukzQSSbDY2NgAoMF7QtSRSrRYJNLS0jB37lx07twZ33zzDTQ0NPDtt99y8w0NDbFo0SJMmjQJ1dXV0NRUqfKJArw8xgJQsBCijlSmxRIZGQl/f3/4+PggODgYOjo6EAqF2LhxI54/f84tV1VVBU1NTfD5fCVWSxqLJFgsLCygpaWF8vJyVFRUKLkqQkh9qESwPH78GLNmzcKCBQuwZMkS8Hg8AICenh5u3LiBLVu2oLy8HDk5Odi8eTPGjh3LLUOaF0mwmJiYwNDQEAC1WghRNyoRLD///DNKS0uxZcsWqetVNm3ahO3bt0MgEGDgwIEYNWoUunbtis8++0zZJZNG8nKwSE7QoNu6EKJeVGKQYsWKFVixYkWd83ft2tWE1RBlqaqqQklJCXg8HoyMjKjFQoiaUokWCyHAvy0TY2NjaGhocC0WOuWYEPVCwUJUhuTiSBMTEwDggoVaLISoFwoWojIk4yumpqYAQF1hhKgpChaiMl4euAf+DRYavCdEvVCwEJXx32AxNjYGQMFCiLqhYCEqg1oshDQPFCxEZfw3WGjwnhD1RMFCVEZdwUItFkLUCwULURnUFUZI80DBQlSGJFgkg/bUFUaIeqJgISqjrhYLBQsh6oWChagMyZX3kgskJS0XuqULIeqFgoWoDEmA1NZiYYwprS5CSP1QsBCVIBaLuWCRtFS0tbWhq6uL6upqlJeXK7M8Qkg9ULAQlVBcXAyxWAwDAwNoaWlx0yWtFuoOI0R9ULAQlfDfgXsJOjOMEPVDwUJUgiQ4JC0UCbpIkhD1Q8FCVEJpaSkAQF9fX2o6nXJMiPqhYCEqQRIsBgYGUtPpKZKEqB8KFqISJMHSokULqek0xkKI+tGUd0WxWIzU1FTk5eWBz+ejdevWaN++PXg8niLrI28I6gojpPmod7BERkbiwIEDuHHjBvdhAAA8Hg9GRkYYOHAgJk2aBFdXV4UWSpq31wULdYURoj5kDpb09HR8/vnnyMrKgo+PD7Zv345OnTrBxMQEjDEIBAI8fPgQkZGRWLx4MWxsbLBmzRrY2to2Zv2kmagrWCQXS1KLhRD1IXOwfPrpp5g3bx6GDBlS63xLS0tYWlrC09MTn376KS5duoRPP/0Ux48fl2n7ERER2Lx5M9LT09GqVSvMmDEDEyZMQFFREVavXo2IiAi0aNECn3zyCd555x1ZyyZqoq5gMTMzAwA8ffq0yWsihMhH5mA5evSozOMnPB4PQ4cOhY+Pj0zLZ2VlYf78+diwYQO8vb0RGxuLmTNnwtraGmFhYdDQ0MD169eRnp6OGTNmwMbGBu7u7rKWTtRAWVkZgJrB0rlzZwBAUlJSk9dECJGPzGeFvRwqK1euREJCQo1lBAIBvL29a13nVZ48eQJfX18MHToUGhoacHZ2hru7O27evInw8HAsXLgQenp6cHBwwLvvvotff/1V1rKJmqirxdKxY0doaGjg0aNHqKioUEZphJB6kut049DQUEyYMAFnzpyRmi4Wi+XqsnBzc8OaNWu49wUFBYiKikKbNm3A4/HQoUMHbl7Hjh2RmJgoT9lEhZWUlACoebqxjo4O2rVrB7FYjLS0NGWURgipJ7mvY/nkk08QGBiIoKAgiMVihRVUXFyMOXPmoEePHnB0dIS2trZUy0dPT4/udNsM1dUVBgBdunQBACQnJzdpTYQQ+cgdLL6+vjh8+DAuXryI6dOnQyAQNPgalrS0NPj7+6N169YIDg6Gvr4+KisrpZ7FIRQKa3yrJeqvrq4wgMZZCFE3cgWLJEAcHBy4s77GjRvXoC6qyMhI+Pv7w8fHB8HBwdDR0UH79u3BGENmZia3XGpqKvdBQ5qPum7pAlCLhRB1I1ewvNyCaNmyJfbt2wcfHx/MmjVLriIeP36MWbNmYcGCBViyZAkXXPr6+vDx8cHmzZtRWlqKhw8f4tixY/Dz85NrP0R11XVLF4BaLISoG7lu6fLxxx9LfQDw+XysXr0ajo6OCA0Nrff2fv75Z5SWlmLLli3YsmULN33SpElYs2YNvvzyS3h5eUFbWxuzZ8/G4MGD5SmbqLBXdYVJWiwpKSkQi8XQ0KBb3BGiynhMxoeJ16cbQpW6qjIzM+Ht7Y3Lly+jbdu2yi6H1KFLly4oKyvDw4cPazyTBQBcXV2RnZ2Nmzdvol27dkqokBAiK5lbLL6+vuDxeGCMSf0XkO4a4/F4ePDggeIrJc2WWCzmzgqr68SMTp06ITs7G8nJyRQshKg4mYPl8uXLUu8ZYxg1ahT27NkDKysrhRdG3hxCoRDAi1PJ+Xx+rct06dIFN27cQFJSEry8vJqyPEJIPckcLNbW1rVOt7CwqHMeIbJ41fiKBJ0ZRoj6oFFQonSyBAudGUaI+qBgIUr3qlONJV4OFhnPNyGEKAkFC1E6WVosFhYWMDAwQEFBAQQCQVOVRgiRg8xjLD///HONaSKRCKdOnYKpqanU9MmTJze8MvLGeNVV9xI8Hg9dunTB33//jaSkJLRq1aqpyiOE1JPMwRISElJjWuvWrWtcEMnj8ShYSL3I0hUGvOgOkwRL3759m6I0QogcZA6WP/74ozHrIG+wV93Z+GWS61eysrIavSZCiPxkDpZhw4bBw8MDHh4e6Nu37yu7LQipD8mzWF4XLJIWjSSICCGqSebB+927d6Nz5844ceIEvLy8MGHCBGzfvh337t1T6PNYyJtHlsF7gIKFEHUhc4vF1tYWtra2mDx5MkQiEf7++2/cuHEDGzZsQHp6Otzc3LgWTfv27RuzZtLMULAQ0rzIdXdjPp8PNzc3uLm5YcGCBSguLsbNmzcRERGB/fv348KFC4qukzRjso6xSOZTsBCi2uQKlv8yNDTEsGHDMGzYMEVsjrxhqMVCSPMi0xhLamoqbt68yQ2ySly5cqVRiiJvFllPN5bMlyxPCFFNrw2Wn3/+GXPnzsWBAwfg6+uLixcvcvOCg4MbtTjyZpC1xaKnpwfg37shE0JU02u7wo4cOYLQ0FC0aNECGRkZWLhwITIzM/HBBx/QPZuIQshy5T1AYyyEqIvXBotIJOK6IGxsbHDgwAEsXLgQz549o2AhCkFjLIQ0L6/tCmvdurXUEyENDAywe/du5OfnIzExsVGLI2+G+o6xULAQotpe22LZuHFjjaf6aWpqYuPGjXjvvfcarTDy5pD1dOOXg+XlR2MTQlTLa1ssFhYWMDMzAwCEhYVJzevVq1fjVEXeKLLe0oXP50NXVxeMMZSXlzdFaYQQOdTreSyBgYGYPHkydYERhRGLxVyL5XVdYcC/Z4ZRdxghqqtewRIaGgoNDQ2MGzcO69evp+sJSINJTh3W09Or0eVaGxpnIUT11StY7OzscPDgQWzYsAHnz5/H22+/jd9//12hBcXExKBfv37c+8rKSnTv3h0uLi7cKyAgQKH7JMoj6xlhEnSRJCGqT65buowcORJeXl7YuXMnPvvsMxw7dgyff/45bG1t5S6EMYbjx49jw4YNUtMTEhJgbGyMiIgIubdNVFd9g4WuZSFE9cn9zHs9PT0sWbIEv/32GzQ1NeHn54dvv/1W7kHV4OBg/PLLL5gzZ47U9Li4ONjb28tbJlFxsp5qLEFjLISoPrmDRaJDhw7Yu3cvtm7dijNnzmDEiBFybWfChAkIDQ1F9+7dpabHx8dDIBBg1KhR6N+/PxYsWIDs7OyGlk1UhLxdYRQshKguue9unJGRgYSEBCQmJnL/zcrKkvtqfHNz81qn6+npwdXVFfPmzYOmpibWrl2LefPm4fjx4/KWTlSIrLdzkaBgIUT11StYPv/8cyQmJiIxMZE7m6dt27ZwcHCAr68v7O3tFd5ttWLFCqn3y5cvR79+/ZCVlQVLS0uF7os0vfp2hdEYCyGqr17B8vDhQ9jb28PPzw/29vbo2rWrzF0Y8tq6dSt8fX3RqVMnAEBVVRUAQEdHp1H3S5qGrFfdS1CLhRDVV69gOXr0aGPVUaeEhATExsZi8+bNAIB169ZhyJAhaNmyZZPXQhSPxlgIaX7kGrw/efIkkpKSIBKJFF1PDevWrYORkRGGDh0KLy8vaGlpYePGjY2+X9I0ZL2di4TkrDC6joUQ1SXX4P2mTZsgEAigpaWFLl26wMHBgXvZ29vL3F9emz59+iAqKop7b2pqyrVWSCaabLAAACAASURBVPMj73Us9LAvQlSXXMESERGB7OxsxMXFIT4+HtHR0Th16hSqq6vB4/EQHx+v6DpJM0VjLIQ0P3Kfbmxubg5zc3N4eXkBALKzs7F48WL4+fkprDjS/NEYCyHNT4MvkJQwNzfHunXrcPXqVUVtkrwB6nu6Md0rjBDVJ1ewCASCWqd36NCBusFIvVCLhZDmR66usP79+8Pc3BwODg7o1q0bunXrho4dOyImJoa7zoQQWch75T0N3hOiuuQKlqNHjyI+Ph7x8fG4du0afvjhB1RWVoLP5+Ozzz5TdI2kGaMWCyHNj1zB4uzsDGdnZ+69SCTCkydPYGxsDGNjY4UVR5o/ecdYKFgIUV0yj7EUFRXVOY/P56Ndu3Y1QqWwsFD+ysgbob6nG0uWo8F7QlSXzMEyZcoU7NmzB8XFxa9dViAQYOfOnZgyZUqDiiPNH3WFEdL8yNwV9ssvv+C7777DoEGD4OrqioEDB6JLly4wNTUFYwwCgQAPHz5EZGQkoqKiMHbsWPzyyy+NWTtRc2KxuN5dYbq6ugCA8vJyiMViaGgo7Ix5QoiCyBws+vr6WLVqFT766CP8+uuvOHPmDB48eMDdL0xTUxPdunXD4MGD8dVXX9X5fBVCJCTdqy1atACfz5dpHQ0NDejp6UEoFEIoFDb63bUJIfVX78F7MzMzfPzxx/j4448hFotRUFAAHo8HU1PTxqiPNGOJiYkAgI4dO9ZrPX19fQiFQpSWllKwEKKC5L6lC/Di2yPdvp7IS3Ixbbdu3eq1Ho2zEKLa5AoWe3t78Hi8WudpaWnB3NwcI0eOxPz582Xu4iBvHgoWQponuYJlzZo12Lp1K+bPn4+ePXsCAO7fv4/g4GD4+/ujS5cu2L59OxhjWLRokUILJs2HvMEieSYLBQshqkmuYAkJCeGe5Chhb28PCwsLrFu3DufPn4eZmRkWL15MwUJqJRKJ8ODBAwD1DxZ67j0hqk2uczVzcnJgY2NTY7qFhQWysrK4/6cLJEld0tLSUF5eDktLy3qf+EFdYYSoNrmCxdXVFZs2bUJBQQE3LT8/H5s3b4arqysA4MKFC+jQoYNCiiTNj7zdYAAFCyGqTq6usLVr1+LDDz/EoEGDYG1tDbFYjKysLHTq1AnBwcG4evUqvv32W2zbtk3R9ZJmgoKFkOZLrmCxtLTEqVOncOvWLSQkJEBTUxNdunRBv379ALzoA7969SqdikzqpIhgofuFEaKa5L6Ohc/nQ0NDA5qamhCLxaiqqkJ1dTU0NTUpUMhrUYuFkOZLrmDJzc3FvHnzEB8fD2trazDG8PTpU9ja2mL//v1o1aqVouskzUh+fj6ysrKgq6sLW1vbeq9PD/siRLXJNXi/bt068Pl8XL58GeHh4bhw4QIuX74MExMTBAUFKbpG0sxIWisODg5yXUBLLRZCVJtcwfLXX39h9erVUjeaNDc3x2effYY///yzQQXFxMRwYzUAUFlZicDAQLi7u6Nv377YvXt3g7ZPlK8h3WAAPZOFEFUnV1eYjo5Orbd04fF43N2O64sxhuPHj2PDhg1S07dt24a0tDRcvHgRxcXFmDlzJszNzTFmzBi59kOUr6HBQlfeE6La5GqxeHh4ICgoCM+fP+emPX/+HOvXr8eAAQPkKiQ4OBi//PIL5syZIzU9LCwMs2fPhrGxMdq2bYsZM2bgyJEjcu2DqAZJsDg6Osq1PnWFEaLa5GqxLFu2DO+//z48PT1hZWUFAHj69Cm6du2KVatWyVXIhAkTsHDhQty+fZubVlRUhNzcXHTu3JmbZmtry91unagfgUCABw8egM/nw8HBQa5t0OA9IapNrmAxMzPD6dOn8ddffyE5ORm6urro1KkT+vfvL3chtT0YTPKNVPLUQOBFN0h5ebnc+yHKde7cOYhEIgwZMgQGBgZybYNaLISoNpmDJTk5ucY0GxsbqXuGSZZ5uYXREJK+9IqKCm6aUCiU+TG2RPWcPn0aAODn5yf3NmjwnhDVJnOw+Pr6gsfjgTHGTfvvAD5jDDwej7trbUMZGxvDzMwMqampXIsmLS1NYcFFmlZubi5u3LgBLS0tDB8+XO7tUIuFENUmc7Bcvny5Meuok5+fH3bs2IGuXbuirKwMISEhmDZtmlJqIQ1z9uxZiMVieHp6wsTERO7tULAQotpkDhZra+vGrKNOCxcuxPr16+Hr6wuxWIz33nsPEydOVEotpGF+++03AA3rBgP+7SKlwXtCVBOPvdy31QxlZmbC29sbly9fRtu2bZVdzhvr2bNncHNzg7a2NqKjo2FoaCj3thhjaNeuHcRiMR49egRNTblveUcIaQRyXcdCSH2dPXsWjDF4eno2KFSAF2N71B1GiOqiYCFN4ty5cwAa3g0mQcFCiOqqdx+CQCDA1atXkZSUhNLSUhgYGMDBwQGDBw+GkZFRY9RImoG0tDQAQK9evRSyPQoWQlRXvYJl37592Lp1K3g8Htq2bQtDQ0OUlJTg4MGD0NTUxKJFi+iMLVKDWCzmbv9jZmamkG1SsBCiumQOlhMnTmDr1q1YtmwZxo0bJ3U1fEVFBUJDQ7Fx40aYm5s36BoF0vzk5+ejuroaJiYm0NHRUcg2KVgIUV0yB8vBgwexdOlSTJo0qcY8HR0dTJw4EWVlZThw4AAFC5GSk5MDQHGtFYCChRBVJvPgfXp6+mvvXOzl5YWUlJQGF0Wal9zcXAAULIS8KWQOlvLy8tfeNNDQ0BCFhYUNLoo0L5IWS203GpWXJFhKSkoUtk1CiGLU63Tj2h7uRcjrNEaLpV27dgCgsPvSEUIUp15nhYWGhr7yzsJ0t1lSG0mLpU2bNgrbZp8+fQAAt27dUtg2CSGKIXOwWFlZ4ZdffnntcpaWlg0qiDQ/jdFicXNzg6amJmJjY1FUVETXUBGiQmQOlj/++KMx6yDNWHZ2NgDFtlhatGiBHj164O7du4iMjIS3t7fCtk0IaRi6pQtpdJIWiyKDBQD69esHgLrDCFE19QqWkpISnDx5kjsThzGGvXv3Ytq0afjkk08QExPTKEUS9dZYwSIZZ7l586ZCt0sIaRiZg+Xx48d4++23sWbNGuTn5wMAgoKCsGXLFrRu3RpGRkZ4//33ce/evUYrlqifiooKFBQUQFNTs0EP96pN7969oaGhgfv379OJI4SoEJmDZevWrXB0dMSNGzdgY2OD/Px8HD58GG+99Ra2bNmCNWvW4OOPP8b27dsbs16iZiT3CGvdujU0NBTb82poaAgnJydUV1fj7t27Ct02IUR+Mv9Lv3nzJubNm8fdI+z69esQiUQYM2YMt4yHhweio6MVXyVRW41xqvHL+vbtC4C6wwhRJTIHS3FxMVq1asW9v3XrFjQ1NeHu7s5N09fXh1gsVmyFRK01xn3CXiYJFhrAJ0R1yBws1tbWSE5OBgCIRCJcu3YNvXv35p4/Drz4x02P/yUva4zbubzM3d0dPB4P//zzD4RCYaPsgxBSPzIHy7hx47B27VqcPn0agYGBeP78OSZOnMjNj4qKwnfffYeRI0c2SqFEPTXGxZEvMzExgYODAyorK/HPP/80yj4IIfUjc7DMmDEDPj4+WLduHa5cuYLFixdj6NChAICvvvoKU6ZMQa9evTBjxoxGK5aon8YeYwEAFxcXAMD9+/cbbR+EENnJfOU9n8/HsmXLsGzZshrz/P39MX78eDg4OCi0OKL+GrvFAgCOjo4AgLi4uEbbByFEdjIHi1AoxNdff43w8HBoa2vDx8cHS5cuhYGBAbp27dqYNRI11hQtlu7duwOgYCFEVcjcFRYcHIwrV65g5syZmD59Oq5du4bVq1c3Zm1Sjh8/DkdHR7i4uHCvsLCwJts/kU9jnxUGAA4ODuDxeEhKSkJ5eXmj7YcQIhuZWyzh4eH45ptvuNM73d3dMXnyZFRVVUFLS6vRCpSIj4/HBx98gKVLlzb6vohiMMYa7XYuL2vRogU6deqE5ORkJCYmwtnZudH2RQh5PZlbLNnZ2ejYsSP33snJCYwx5OXlNUph/xUXF0djOGqmqKgIFRUVMDAweOVzfBSBusMIUR0yB4tIJAKfz+fe83g8aGlpoaqqqlEK++++ExIScOrUKQwYMABDhw7Fnj17wBhr9H0T+TXFwL2EZAA/Nja20fdFCHm1ej1BUlkEAgG6d++OMWPGYPv27UhJScHcuXOhr6+PyZMnK7s8UoemGLiXoBYLIaqjQY8mFolEOHXqFExNTaWWU/SHvZmZGQ4dOsS9d3BwwJQpU3DhwgUKFhWmjBZLfHw8xGKxwm94SQiRXYMeTdy6dWuEhoZKTePxeAr/sE9KSsK5c+ewYMECblpVVRV0dHQUuh+iWI19O5eXtWrVChYWFnj27BnS09OlxgMJIU1LLR5NbGRkhH379sHCwgLjx49HfHw8Dh48iMDAQKXVRF6vKU41fpmjoyOePXuG2NhYChZClEgt+gvMzc2xc+dOHDlyBL169cKCBQswd+5cvPXWW8oujbxCU46xADTOQoiqUIvBe+DF883/2+1GVFdeXh6uXLkCALCxsWmSfdKtXQhRDWrRYiHqZ+XKlcjLy0P//v25i2obG7VYCFENFCxE4U6fPo0zZ85AX18fW7ZsabIztGxsbGBoaIicnBw8efKkSfZJCKmJgoUoVG5uLlauXAkAWL16dZN1gwGAhoYGBg4cCABSp6cTQpoWBQtRqG3btiE/Px8DBw7E1KlTm3z/H330EQDgwIEDKC0tbfL9E0IoWIiCSU5LX7JkCXg8XpPvv3fv3ujVqxcKCgpw5MiRJt8/IYSChShQRkYG0tLSYGRkxD3VURnmzJkDANi7dy+qq6uVVgchbyoKFqIwf/75JwDAw8MDmprKO5N92LBhsLW1RUZGBs6ePau0Ogh5U1GwEIW5du0aAGDQoEFKrYPP53NjLd9//z3dBZuQJkbBQhRCJBIhIiICgPKDBQDeffddmJqa4v79+3RdCyFNjIKFKERMTAwKCgrQvn17dOjQQdnlQE9PD35+fgBAj7AmpIlRsBCFUJVusJeNHTsWAHDy5EmIxWIlV0PIm4OChciMMYZjx47hzp07Neb99ddfAFQrWHr16oW2bdvi2bNnuH37trLLIeSNQcFCZBYeHo5PPvkEY8eOxQcffIDk5GQAQElJCaKioqChoQEPDw8lV/kvDQ0NjBkzBgB1hxHSlChYiMxOnDjB/f+FCxfg5eUFDw8PjBo1CtXV1ejZsyeMjY2VWGFNkmA5e/YsKisrlVwNIW8GChYik8LCQly+fBk8Hg/nz5/HlClTAADp6elITEwEAJV8Po6DgwPs7e1RUFCAq1evKrscQt4IFCxEJufOnUNFRQX69+8PJycnbNiwAXFxcfjzzz9x6tQpHD9+nLt2RNVIWi3Hjx9XciWEvBkoWIhMJA9ZGzduHDfN0NAQnTp1gpubG/r16wctLS1llfdKY8aMgYaGBs6ePYslS5ZAKBQquyRCmjUKFvJaWVlZuHHjBrS1tfH2228ru5x6s7GxwebNm6Grq4sjR47A19cXf/zxB3JycuiqfEIaAQULea3Tp0+DMQYfHx+VG5yXlb+/P3777Td07NgRDx8+xNSpU+Hi4oKePXtKnZRACGk4ChYAVVVV9M31FSSn6krGKtRVt27dcP78eXz88cdwd3eHkZERnj9/jmXLliEzM1PZ5RHSbLwxwSIUCjF69Gh4enoiKCgIUVFROHbsGCZNmoSOHTviiy++UHaJKik5ORn379+HoaEhvL29lV1Og+nr62PFihUICwtDfHw8Ro8ejfLycvzvf/9TdmmENBs81sy/qmdmZsLb2xvDhg3jBqBro6GhgStXrqBz585NWJ3q27x5M7Zs2YL33nsPW7ZsUXY5CpeVlYXBgwejtLQUhw4dgqenp7JLIkTtvTEtlrCwMGhpaeG7775DQEAAbG1t0adPHwQFBeGdd96BWCzG5s2ba12XMYbCwsImrlj5GGM4efIkAGD06NFKrqZxWFpaYtGiRQCA1atXo6KiQskVEdIMMDXx4MED5u/vz3r06MF8fX1ZdHS0TOtlZGQwOzs7Zm1tzbZu3VrrMpmZmaxDhw7MysqKxcXF1Zj/9ddfMysrKxYQEMCSkpIadBzqJCYmhllZWTFnZ2dWVVWl7HIaTUVFBRs8eDCzsrJib731FgsJCWG5ubk1lsvPz2dpaWlNXyAhakYtWiyVlZWYO3cu3n77bURGRmL27NmYMWMGSkpKZN6Go6Mj5s6dW+s8a2trTJ06FQCwadMmqXnZ2dnYu3cvAOD8+fPw9PTEihUr3ohvtqdOnQIA+Pr6KvWJkI1NW1sbmzZtgpGREWJiYhAYGAhXV1eMHz8eu3fvxsWLFzFv3jy4urpiwIABNZ5KKRAIkJubW+f2c3NzsWLFCly/fr2xD4UQlaAWYyzXr1/HihUruDvoAsDEiRMxduxY+Pv7v3JdyRhLSEgIBgwYUOdyubm56NevH4RCIX777Te4uroCAL788kvs2bMHXl5esLa2xuHDhyESibBo0SIsXbpUMQeogsRiMfr06YOnT5/i5MmT6N27t7JLanRCoRAXLlxAaGgorl69iurq6lqX09fXx5kzZ2BnZ4fIyEhMnz4dRUVFGDJkCCZNmgQfHx/uYlGhUIh33nkH0dHR0NXVRVhYGJydnZvysJSCMYYHDx7g9u3buHXrFp49e4aWLVuidevWaN26NSwtLWFhYYGOHTvC1tYWfD6/UeqorKzEP//8g+LiYhgZGcHIyAiGhoYwNDSEvr4+NDTq/m4tFotRVlaG4uJiFBcXo6CgANnZ2cjNzcXz589RUFCAgoIC8Pl8WFlZwcrKCjY2NmjXrh1sbGzA5/NRUlKC0tJSmJiYQF9fv9615+bm4unTp3j69CmAF123lpaWEIlEyM3NRU5ODgQCAVdLcXExSktLIRQKwefzoaWlBV1dXZibm8PS0hJt2rSBnp4edHV1UVFRgdTUVKSmpiIzMxO5ubnIzc0Fj8dDy5Yt0apVKxw4cECun7taBMv+/ftx9epV7N+/n5u2atUq6OnpYfXq1a9cVxIsly9fRtu2bV+5bFBQELZv344OHTrg9OnTAIA+ffpAKBQiPDwc3bt3R0REBPz9/aGlpYXw8HB07dq1xnYyMjKQkJCA4uJilJWVYdCgQbCxsan/gSvRnTt3MHbsWFhZWeH27duv/AfYHBUVFeHKlSu4ePEiUlNT4ePjA39/f3z99dc4deoUbG1tsXjxYnz66acoLy+XWrd9+/b48ssv4eXlhVmzZuHcuXPQ0tJCVVUVLCwscPbsWVhYWCjpyBpXbm4ujh49isOHDyM9PV2mdfT09ODo6Ag3NzcMHjwY7u7u0NXVBfAioNLS0nDt2jVER0fDzMwM7dq1Q7t27WBhYQFzc3Noa2sjJSUFSUlJyMnJQXl5OSoqKvDgwQPcuHEDpaWlde5bW1ube0n+xiWB8t/fa0MZGBigTZs20NbWBo/Hg4aGBnR1ddGiRQtoaWlBKBRCKBSiuLgYeXl5KCgoUOj+5fHkyRO51lOL/o2ysjLuD01CT09P4bfmWLBgAa5du4b79+/j/fffh6urK4RCIXx8fNC9e3cAgIeHByZPnoyff/4Zy5YtQ1hYmNQf5M6dO7Fp0yapb7tGRkY4dOgQevXqpdB6FW3z5s24d+8e+vbti/v37wN4MWj/poUK8OJ3Nnr06BonLXzzzTdITEzEgwcPMH/+fADA5MmTsXTpUoSFheHgwYNIS0vD9OnT0alTJ6SkpMDY2BjHjh1DYGAgbt++jYCAAKxZswZ6enpgjCEhIQFxcXEQCATw8fGBj49Pjb/3piQUCnH37l3cuXMHUVFRqKys5L7tm5iYcN++k5OTERsbi6SkJIhEIvB4PFRUVHAPVTM3N4eHhwf69u2Ljh07Ij8/H8+fP0dubi6ysrKQlZWFxMREPH36FFFRUYiKisKuXbugq6sLU1NTaGpqorKyEtnZ2Q06Hjs7O1hZWaGoqAhFRUUoKSnhXpWVla+867Wuri6MjY1haGgIY2NjtGnTBm3atEHr1q1hamoKY2NjVFVVca2KjIwMPH78GJmZmWCMwdDQEHp6eigoKOD2KSsNDQ2YmZnBysoKlpaWAMD93LS0tGBmZgYzMzO0atWK+70YGBjAwMAALVq0gFgsRlVVFUpLS5GdnY2srCzk5uaivLwc5eXl0NTUhK2tLTp27AgbGxvuuHg8HvLy8pCfny/3z1xtWizXrl3Dvn37uGmrVq1CixYtsGrVqleuW58WC/BiTMXPz0/qgrmXu8aAF3f69fT0RHZ2Nr744guMHDkSBQUF+PLLL7nnvnt4eKBVq1bIzs7G7du3oa+vj4MHD6JPnz71Pfwm8eeff2LixIk1pktaauRf6enpGDFiBAoLC/HJJ59g6dKl4PF4AF5cbLtv3z5s3rwZJSUl0NTUxOHDh+Hh4YG8vDyMHDkSGRkZr9y+kZERxowZg08//RQtW7ZslGMoLCxEdHQ0CgsL4erqCmtrawgEAvz444/48ccf5T4Lks/nw8fHB5MmTYKnp6dMXVwCgQAxMTGIiIjA1atXER8fLzXfxMQEgwYNQp8+fVBYWIhHjx4hIyMD2dnZyM7ORnl5OWxtbdGlSxe0bdsWurq60NbWhpWVFQYOHAgrK6ta9ysWi7lgqays5C6S5vF40NPTg56entxfqkQiETQ0NLi/C8YYCgoKkJubi+rqajDGIBaLIRQKUVZWhsrKSujp6aFFixYwNDRE69atYWJiorZf6tQiWP766y8EBgZK3fZ84sSJGDNmDN57771XrlvfYAGAxMREjBkzBoWFhRg0aBB++eWXGsucPXu21rv5tm7dGt999x13PUR1dTU++eQThIWFQU9PDwcPHkS/fv1kqqOxlJSU4MmTJ1w3XkVFBXx8fJCamooJEyaAz+fjxo0bcHR0xK5du7h/HORfGRkZePbsWZ1jTzk5OQgJCUGfPn3g5eXFTU9OTkZQUBD3gSgWi9GpUyc4OjpCR0cHp06d4lqLrVu3xoYNG+p8HIFIJEJERAQuXbqE9u3b47333oOBgUGN5fLz8xEdHY24uDjExcXh/v37SE1NlVrG2toa+fn5KCsrA/DicQP9+/dHnz59YGJigqKiIhQWFqKgoAD5+fkoLi5Gu3bt4OzsDHt7e671xefzoaOjI9fPVELy7V4kEgEA2rZt+8qAEovFavsB3FypRbBUVlZi6NCh+OCDDzB58mRcuHABgYGBuHTp0mu/0ckTLAAQGRmJ7du3Y+XKlbWOozDGsHz5cvz+++/Q0dGBrq4unJyc8OWXX6JNmzZSy4pEIixduhRHjx6Fubk5/vjjD5iYmMhci6KkpKRg//79OHr0KEpKSjBp0iR89dVX2Lt3L9avX49OnTrh0qVL0NbWbvLayL8ePnyI1atX4+bNmwBenJU3d+5c9OjRA8CLcDp27BiOHz+OZ8+ecesZGhpiwoQJsLS0REFBAXJycnDv3j3ueTkv09HRgaOjI4yNjXH37l0UFRUBAIYMGYIFCxaobMuaqAe1CBbgRSviiy++wMOHD9G2bVusXLlSpm/+8gaLoolEIowbNw5RUVF499138d133zXp/kNCQvD5559z7/l8PkQiERwdHZGSkoLy8nIcOXIEAwcObNK6SO3EYjH279+PdevWcYPIPXr0gLa2NiIjI7nl2rdvj1GjRiEqKgq3bt2qdVs6OjpwdnaGk5MTHB0d0a1bN9jb23NfIMRiMRISEqClpUV3niAKoTbBIi9VCRbgxTfN4cOHo7y8HD/99BN8fHy4eWKxGKmpqRAIBOjdu7dCu58iIyPxzjvvQCQSYcKECQgICABjDB999BEePXoEAPDz88P333+vsH0SxcjIyMC+ffvw66+/cmcJ6evrY9SoUfD394e7uzv3txIdHY3Q0FAwxmBiYgJTU1M4OTnBycmpwd1ThNQHBUsT27VrF7766itYWFhg7ty5SElJQUJCAmJjY7kzRiRnDSkiXPLz8zFs2DA8ffoUs2fPRmBgIDevsLAQK1euRFJSEn766SfuzBOieiSnvDPGMGzYsHpfE0FIU6JgaWIikQhjxozBvXv3asyztLREXl4eKisrMW3aNKxbtw4AcO/ePQgEAjg7O9fr+gfGGAICAnDhwgW4uLggNDSUxk8IIY1OLa5jaU74fD62b9+OtWvXwtTUFJ06dUKXLl3QvXt3tGnTBpcvX8aHH36IAwcOICEhAampqVK3C7GwsMDIkSOxatWq13ZvHDx4EBcuXICRkRG+//57ChVCSJOgFosKunr1KmbMmMEN2trY2MDGxgb3799HcXExAMDd3R0hISF1nhWXlZUFT09PFBcX4/vvv4efn1+T1U8IebNRi0UFDRkyBMeOHcOtW7cwePBgdOvWDTweD2KxGPfu3cOsWbNw584d+Pn54ZtvvkGnTp24K2YlAgMDUVxcjOHDh2PUqFFKPBpCyJuGWixqKCsrC++//z7i4uK4abq6uhg8eDDmzJmD58+fY+bMmdDX18fVq1frvPKYEEIaA7VY1JClpSVCQ0Oxbt063L17F0+ePEFBQQHCw8MRHh7O3Vl3xYoVFCqEkCZHwaKmDAwMEBQUxL3PycnBvn378NNPP6GwsBAuLi6YNm2aEiskhLypqCusmSkpKcGVK1fg4eHRaDcwJISQV6EWSzNjYGBAg/WEEKWiW4ISQghRKAoWQgghCkXBQgghRKEoWAghhCgUBQshhBCFomAhhBCiUM3+dGPJc7NffoQrIYQQ2VhYWEBTs35R0eyDRXLL+cmTJyu5EkIIUT/yXFze7K+8Ly8vR2xsLMzMzMDn85VdDiGEqBV5WizNPlgIIYQ0LRq8J4QQolAULIQQQhSKgoUQQohCUbAQQghRKAoWQgghCkXBQgghRKEoWAghhCgUBYsKOn78OBwdHeHiXB2XfAAAC/VJREFU4sK9wsLClF2WwsTExKBfv37c+8rKSgQGBsLd3R19+/bF7t27lVid4tR2nN27d5f6vQYEBCixQvlFRERg3LhxcHV1xdChQ3HkyBEAQFFRERYsWIBevXph4MCBOHHihJIrbZi6jjMjIwP29vZSv8tVq1YpuVr5XblyBaNGjYKLiwt8fHy445T73yYjKufLL79kmzZtUnYZCicWi9nRo0dZr169WK9evbjp33zzDZs8eTIrKChgGRkZbPjw4SwsLEyJlTZMXccZExPD+vfvr8TKFOPp06fMxcWFXbhwgYlEIhYdHc169+7N/vzzT7Zo0SK2cOFCVlZWxuLj41m/fv3Y7du3lV2yXF51nOfOnWPvvvuusktUiOzsbNa9e3d29epVxhhjsbGxzMnJicXGxsr9b5NaLCooLi4ODg4Oyi5D4YKDg/HLL79gzpw5UtPDwsIwe/ZsGBsbo23btpgxYwb3jUkd1XWccXFxsLe3V1JVivPkyRP4+vpi6NCh0NDQgLOzM9zd3XHz5k2Eh4dj4cKF0NPTg4ODA9599138+uuvyi5ZLnUd571795rN7xIA2rRpg5s3b2Lw4MEQi8UoKCgAn8+Hvr6+3P82KVhUjEgkQkJCAk6dOoUBAwZg6NCh2LNnD1gzuPPOhAkTEBoaiu7du3PTioqKkJubi86dO3PTbG1tkZiYqIwSFaK24wSA+Ph4CAQCjBo1Cv3798eCBQuQnZ2tpCrl5+bmhjVr1nDvCwoKEBUVhTZt2oDH46FDhw7cvI4dO6rt77Ku4+zWrRvi4+ORlJSE4cOHY8CAAVi1ahWKioqUWG3DGBgYQCgUwsnJCQEBAZg8eTJatmwp979NChYVIxAI0L17d4wZMwZ//PEH9+338OHDyi6twczNzWtMKysrAwDo6upy0/T09FBeXt5kdSlabccJvDguV1dX/PTTTzh//jx0dXUxb968Jq5OsYqLizFnzhz06NEDjo6O0NbWBo/H4+ar++9S4uXj9Pb2hpGRETw8PHD8+HGEhYUhKysLgYGByi6zQXR0dPD333/j+PHjOHHiBA4cOABAvn+bzf62+erGzMwMhw4d4t47ODhgypQpuHDhQrO89b+enh4AoKKigpsmFArRokULZZXUaFasWCH1fvny5ejXrx+ysrJgaWmppKrkl5aWhrlz56Jz58745ptvkJKSgsrKSjDGuHBpDr/L/x6nhoYGvv32W26+oaEhFi1ahEmTJqG6urredwJWFRoaGtDW1oaTkxP8/f0RGxsLQL5/m9RiUTFJSUkIDg6WmlZVVQUdHR0lVdS4jI2NYWZmhtTUVG5aWlqaVPO7udi6dStSUlK491VVVQCglr/byMhI+Pv7w8fHB8HBwdDR0UH79u3BGENmZia3XGpqqlr/Lms7TqFQiI0bN+L58+fcclVVVdDU1FTLR3PcuXMH48aNk5pWWVkJIyMjuf9tUrCoGCMjI+zbtw9Hjx6FWCxGbGwsDh48WOMX35z4+flhx44dEAgEyMzMREhICPz8/JRdlsIlJCRg/fr1KCoqQlFREdatW4chQ4agZcuWyi6tXh4/foxZs2ZhwYIFWLJkCdc60dfXh4+PDzZv3ozS0lI8fPgQx479X3v3G9LU28YB/IvZMM4k51glJlv1YswSWhZKiFOr0R+IiqIMptGWtgyRRLTIfzVqyyW12NIoEq0In5xIVpRSlNZYlNmLhBLqRX9Ai82IxbTNPS+kg0tz+9V5Hn+b1wfOi+3sXPe1I2eX9322+/5PyP4tf/c+58yZgydPnqC2thZutxuDg4M4ffo0tm7d6jcMGCpkMhkGBgZw+fJleL1e9PT0oKWlBdu3b//ja5PWY/kXstlsqKmpwbt37yAQCKBWq8NqGMxut6OgoADPnj0DMNbV1uv1uHv3LkZHR7Fz504UFRWF5EU63q/v0+l0QqfTobu7G16vFwqFAhUVFZg7d+40Z/rPnDx5Eg0NDROGRHbv3o28vDxUV1fj8ePH4PF42LdvH3JycqYp078z1fvMzs6GTqfD8+fPERERgU2bNqG0tDQke5/A2DcWdTod3rx5g7i4OBQWFkKpVP7xtUmFhRBCCKdoKIwQQginqLAQQgjhFBUWQgghnKLCQgghhFNUWAghhHCKCgshhBBOUWEhIS8rKwtSqXTSzePxTHd608put7Pn4tixY/jw4QOkUumkEwlarVakpKQEjGk0GtmYDx48+F+kTUJcaE5qQ8gvDh06NOnsBKE6bxPXOjs7IRAIMDQ09NextFotcnNzkZaWxkFmJBzRVUfCAsMwEIlE053Gv1ZsbCwYhuGksDAMA4ZhOMiKhCsaCiNh79y5c9BoNFCr1UhOTmaXeb506RIyMzMhl8uRnZ2N3t5e9hiPxwODwYDU1FSkpKTgwoUL2LZtG6xWKwBApVLBYDD4tfPr0NBU8cvKylBVVYXDhw9DLpcjKysLFouF3e/1emE2m5GRkQG5XA6VSoX+/n4MDg4iMTERNpvNr+1169Zxugyw1WqddGhRpVJx1gYJX1RYyIzQ1dWFlStXorm5GQqFAtevX0djYyMqKyvR2toKhUKB3NxcvH//HsBYMWpra4PBYEBDQwPu37+P/v7+oNsLFB8Abty4gfnz57MT/p09e5adqtxsNqOpqQlHjhxBa2sr4uLikJ+fD6FQiNTUVLS3t7Nxent7MTAwAKVSydHZAjZu3Iju7m52O3PmDCIjI6HRaDhrg4QvKiwkLOj1esjlcr/Nbrez+6OiopCfn48lS5YgNjYW9fX1KC4uRkZGBiQSCfbv34/k5GRcu3YNPp8Pzc3N0Gq1UCgUkMlkMBqN8Hq9QeczVfyfEhISUFRUhMWLF+PAgQOIiYnBq1ev4PP52KWNlUolJBIJKioqoFQq8fXrV2zevBkdHR3stPvt7e3IyMhAdHR00Pnt2LFjwvmqrKz0O18ikQgikQgejwfHjx9HQUEBFApF0G2QmYvusZCwkJ+fP2E67/ErOS5cuBAREWP/R7lcLnz69Anl5eV+H6YjIyPg8XhwOBxwOBxYunSp3/ELFiwIKpdA8X8Si8V+xzEMA4/HA6fTCYfDgaSkJHYfn89HWVkZgLFhr6qqKnR3dyM9PR137txBdXV1ULn9ZDKZ/JYQBoB79+7h4sWLfs8NDw/j4MGDWL58ObRa7T9qg8xcVFhIWBAIBBM+qMcbP5356OgogLFeTmJiot/roqKi/JZiHW98UfjV+K81B4o/VTyfz4fZs2cDwG+nJmcYBmvWrMHt27fB4/EwMjKC9PT03+Y2mbi4uAnnSygUTnhdeXk5XC4XTp06FfLLGJD/HxoKIzNOdHQ0RCIRBgYGIBaL2e3KlSvo6uoCwzCIj4/3u9nudDrx8eNH9jGPx8O3b9/Yx+PvnQSKH0x+QqEQfX197HNutxurV6/Gy5cvAYwtjvbo0SN0dnZi/fr1Uxa9P9XQ0IDOzk6YzWbw+XzO45PwRT0WMiNpNBpYLBbMmzcPSUlJuHnzJq5evYqmpiYAQF5eHmpra5GQkIBFixbBaDRiZGSEPX7ZsmVobGzEli1bwOfzodfr/T7cA8UPZM+ePbBYLIiPj4dEIkFdXR2io6Mhk8kAAGlpaYiMjERLS8uE4Ssu2Gw2GI1GnDhxAjExMfj8+TO7j77WTQKhwkJmpJycHLjdbtTU1ODLly+QSCQwmUxITk4GAOzatQtDQ0OoqqrC8PAw1Go1Xrx4wR6/d+9evH37Fmq1GgKBAIWFhX49mkDxA1Gr1XC5XDh69ChcLhdWrFiB+vp6tnjNmjULGzZsQEdHB1atWsXhmRnT1taGHz9+oKSkZMK+169fc94eCS+0giQhQUpJSUFpaemkv/CfDoWFhRCLxSguLv7ta+x2O3JyctDT08P5jxqlUinq6uqQmZnJaVwS+qjHQkiIefr0Kfr6+vDw4UO/37NMxeFwwOfzcXKvxOVy4fv3738dh4QvunlPSIi5desWTCYTSkpKkJCQENQxa9euRW1tLSftnz9/nuYJI1OioTBCCCGcoh4LIYQQTlFhIYQQwikqLIQQQjhFhYUQQginqLAQQgjhFBUWQgghnPovncm1dN+P2CkAAAAASUVORK5CYII=\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": [
    {
     "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}