ERSP_ITC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import os\n",
    "import json\n",
    "import pyreadstat\n",
    "import datetime\n",
    "import sys\n",
    "import ptitprince as pt\n",
    "import seaborn as sns\n",
    "\n",
    "# configure backend here\n",
    "matplotlib.use('Agg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## File Paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## File paths\n",
    "# CNT import path\n",
    "import_cnt = '[INSERT PATH TO CNT FILES]'\n",
    "\n",
    "#Project name string\n",
    "project_name = '[INSERT PROJECT NAME]'\n",
    "\n",
    "# home directory: where this script is located\n",
    "home_dir = os.getcwd()\n",
    "\n",
    "# Set stdout and stderr to a log file\n",
    "now = datetime.datetime.now()\n",
    "log_file = open(f'{home_dir}/log_{project_name}_{now.strftime(\"%Y-%m-%d_%H:%M\")}.txt', 'a')\n",
    "sys.stdout = log_file\n",
    "sys.stderr = log_file\n",
    "\n",
    "\n",
    "#Make a function to make a directory if it doesnt exist and return the path\n",
    "def make_dir(path):\n",
    "    if not os.path.exists(path):\n",
    "        os.makedirs(path)\n",
    "    return path\n",
    "\n",
    "epochs_path = make_dir(f'{home_dir}/epochs_{project_name}')\n",
    "all_epochs_path = make_dir(f'{epochs_path}/all_epochs_{project_name}')\n",
    "ersp_path = make_dir(f'{home_dir}/ersp_{project_name}')\n",
    "itc_path = make_dir(f'{home_dir}/itc_{project_name}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Epoch event dictionary: maps custom event names to event codes specified in the CNT file\n",
    "epoch_event_dict = {\n",
    "    'con_1': 1,\n",
    "    'con_2': 2,\n",
    "    'con_3': 3,\n",
    "    'con_4': 4,\n",
    "}\n",
    "\n",
    "freqs = [1, 30] # Frequency range for ERSP and ITC values: ex 1-30Hz (1,31) in logspace\n",
    "vmin, vmax = -1, 1.5  # set min and max ERDS values in plot\n",
    "baseline = (-.5, 0)  # baseline interval (in s)\n",
    "wtimes = [0, 1] # time window in seconds \n",
    "morlet_win = [1, 0.5] # number of cycles in the proportion of the window\n",
    "kwargs = dict(\n",
    "    n_permutations=100, step_down_p=0.05, seed=1, buffer_size=None, out_type=\"mask\"\n",
    ")  # for cluster test\n",
    "\n",
    "padding = 1 # when cycles != 1, frequency spacing is divided by padratio (padding)\n",
    "\n",
    "# baseline = [-500 0] # if 'NaN'; no baseline, in some instances this is recommended\n",
    "# morlet_win = [1 0.5] # number of cycles in the proportion of the window;\n",
    "# padding=1\n",
    "# wtimes=[0 1000] # Can be written as [100:100:1000] if you want to find the consequitive windows of 100ms each"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnt_list = os.listdir(import_cnt)\n",
    "cnt_list = [x for x in cnt_list if x.endswith('.cnt')]\n",
    "cnt_list.sort()\n",
    "\n",
    "# Subject list\n",
    "subject_list = [x.split('_')[0] for x in cnt_list]\n",
    "subject_list = list(set(subject_list))\n",
    "subject_list.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Epoch data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Iterate through filelist\n",
    "for subject in subject_list:\n",
    "    # Get file for this subject\n",
    "    file = [x for x in cnt_list if x.startswith(subject)][0]\n",
    "    raw = mne.io.read_raw_cnt(os.path.join(import_cnt, file), preload=True)\n",
    "    \n",
    "    # Set epochs \n",
    "    events = mne.find_events(raw)\n",
    "    reject_criteria = dict(eeg=200e-6)  # 200 µV\n",
    "    flat_criteria = dict(eeg = -200e-6)  # -200 µV\n",
    "    exclude_channels = ['HEOG', 'VEOG', 'M2', 'EKG', 'Startle']\n",
    "    keep_channels = [ch for ch in raw.ch_names if ch not in exclude_channels]\n",
    "    epochs = mne.Epochs(\n",
    "        raw, \n",
    "        events, \n",
    "        event_id=epoch_event_dict, \n",
    "        tmin=wtimes[0]-1, \n",
    "        tmax=wtimes[1]+1,\n",
    "        baseline=(None, 0), \n",
    "        picks=keep_channels,\n",
    "        reject_tmin=-0.5, # -500 ms\n",
    "        reject_tmax=1, # 1000 ms\n",
    "        reject=reject_criteria,\n",
    "        flat=flat_criteria,\n",
    "        reject_by_annotation=True,\n",
    "        preload=True\n",
    "    )\n",
    "    \n",
    "    # Perform baseline correction\n",
    "    epochs.apply_baseline((None, 0))\n",
    "\n",
    "    # Save epochs\n",
    "    epochs.save(f'{all_epochs_path}/{subject}-epo.fif', overwrite=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch_list = os.listdir(all_epochs_path)\n",
    "epoch_list = [x for x in epoch_list if x.endswith('-epo.fif')]\n",
    "epoch_list.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make epochs equal length\n",
    "# Iterate through filelist\n",
    "for subject in subject_list:\n",
    "    # Get file for this subject\n",
    "    file = [x for x in epoch_list if x.startswith(subject)][0]\n",
    "    # Load epochs\n",
    "    epochs = mne.read_epochs(os.path.join(all_epochs_path, file), preload=True)\n",
    "    # Find condition with least epochs\n",
    "    min_epochs = len(epochs[epoch_event_dict[0]])\n",
    "    for condition in epoch_event_dict:\n",
    "        if len(epochs[condition]) < min_epochs:\n",
    "            min_epochs = len(epochs[condition])\n",
    "            min_condition = condition\n",
    "    \n",
    "    for condition in epoch_event_dict:\n",
    "        if condition != min_condition:\n",
    "            current_con_epochs = epochs[condition]\n",
    "            current_epoch_count = len(current_con_epochs)\n",
    "            \n",
    "            # Keep every nth epoch, where n is the number of epochs in the condition with the least epochs\n",
    "            current_con_epochs.drop(np.linspace(0, current_epoch_count-1, min_epochs), reason='Make epochs equal length')\n",
    "            \n",
    "            # Make a directory called {condition}_{project_name} if it doesnt exist\n",
    "            if not os.path.exists(f'{epochs_path}/{condition}_{project_name}'):\n",
    "                os.makedirs(f'{epochs_path}/{condition}_{project_name}')\n",
    "            con_epochs_path = f'{epochs_path}/{condition}_{project_name}'\n",
    "            \n",
    "            # Save separate epochs for each condition\n",
    "            current_con_epochs.save(f'{con_epochs_path}/{subject}_{condition}-epo.fif', overwrite=True)\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform ITC and ERSP Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "morlet_win = [1, 0.5] # number of cycles in the proportion of the window\n",
    "freqs = [3, 30] # Frequency range for ERSP and ITC values: ex 1-30Hz (1,31) in logspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 5.])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tfr_morlet_eeglab_n_cycles(morlet_win, freqs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "30* .7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement tfr_morlet EEGLAB parameter functions\n",
    "def tfr_morlet_eeglab_n_cycles(morlet_win, freqs):\n",
    "    #\n",
    "    if morlet_win[1] < 1:\n",
    "        new_cyc = [morlet_win[0], morlet_win[0]*freqs[-1]/freqs[0]*(1-morlet_win[1])]\n",
    "        new_cyc_list = np.linspace(new_cyc[0], new_cyc[1], int(freqs[1]))\n",
    "        return new_cyc_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.        ,  1.28571429,  1.57142857,  1.85714286,  2.14285714,\n",
       "        2.42857143,  2.71428571,  3.        ,  3.28571429,  3.57142857,\n",
       "        3.85714286,  4.14285714,  4.42857143,  4.71428571,  5.        ,\n",
       "        5.28571429,  5.57142857,  5.85714286,  6.14285714,  6.42857143,\n",
       "        6.71428571,  7.        ,  7.28571429,  7.57142857,  7.85714286,\n",
       "        8.14285714,  8.42857143,  8.71428571,  9.        ,  9.28571429,\n",
       "        9.57142857,  9.85714286, 10.14285714, 10.42857143, 10.71428571,\n",
       "       11.        , 11.28571429, 11.57142857, 11.85714286, 12.14285714,\n",
       "       12.42857143, 12.71428571, 13.        , 13.28571429, 13.57142857,\n",
       "       13.85714286, 14.14285714, 14.42857143, 14.71428571, 15.        ])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(tfr_morlet_eeglab_n_cycles([1, .7], [1,50]))"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABDYAAACsCAYAAACAYSx6AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAAHwBSURBVHhe7d0HfFTF2gbwZ3dTSKUkIaGG3nsJvRcpglhBQUVpYgNFsYECFkQRxYti9179FKUr0ptSQwu9hJ5CSQikkELKbr55Z3chCQkpBEjw+ftbSc5uTpkz5+zMe6YYPDw80kFEREREREREVAwZbf8SERERERERERU7DGwQERERERERUbHFwAYRERERERERFVsMbBARERERERFRscXABhEREREREREVWwxsEBEREREREVGxxcAGERERERERERVbDGwQERERERERUbHFwAYRERERERERFVsMbBARERFRMeeB+v1HYcyYkbi3npttGRUHbnXuxxuzF2BtYBD27NqKdQtnYXhTV9u7RER5Y3J2dp5k+5mIiIiIihUjvFs/gRdHPoj2lRJxZP8ZJKbb3vo3cWyAUV99ibH9OqJW3Cr8uuk8LLa3KC/uTD4y+g/Bd4umY2DTyvB2MyI1zQllqzeAz5k/8Mu2i7ZPERHlji02iIiIiIorYwU88OpkvDBsGJ6fMg79fFi0K9KMpeDfqBma1PSGk21RkXBH8pEjWj71DDp7mWCJWI7xXRugfsP6qN+mP8YvOGX7DBFR3vDbj4iIiKi4ssQjMuoyLEiHJeY8zseznUJR5txjMpYsXYK/Vs3EI162hUXBnchHxspo2bwSTDAjat0vmH8iWS20IDF8Lw6Ep1o/Q0SURwxsEBERERVb0Vj8fAe07dwFbTu9hOUJtsVUNBmNqiIPGAwmGA1FqRh+B/KRyQdlvQzqh3REnY+A2bqUiKhAGNggIiIiulNURffGhTFVEc6ttJYag7BjxxAem/eqoTHX7eZTYa8vO4W6DUnXQtxjtS7HXE/UrSPn84aMptzTLt/5SJ2Pm0lDgxOcbf1xzGlp1h/yRLZr+/EGjOp85H/v8rZuIip6eOkSERER3VJGeLcbiSnTpuLl3pXVb46o2PVFzFqyDUdOheD0sSCs+vEt9K/pYvs84FK9D179Zjl2BZ9CyOlj2L3mv3izXw1cP1eEWnfb4Zj04Yd4Z1hreGUs2Rm90W7kFEyb+jJ6V1ZvmPzQfuR0zNt4AMdOh+JU8C6s/P419K1+bbv54ejXGkPf/RFLAw/ipFpfSMhR7N34J7595wm08pV2CWrfOoxSxz0N74/rgyo5ljrd0Pixifhw2vt4SafPNblvIx/caqHvyzMxd10Qgk+fRkjoKQTvWotfPxqFjhUdbR/KI49auGfUJMyeuwJb9hzGCXUeT4WcxokDW7H0h8kY0swrw3EY4dX+GXzw8XR8MKQJXKSRgrE67nv7I0yfPh3TP/4YU55oms25zSLr+XSsiK4vzsKSbUf0to8FrcKPb/XHtWzkgup9XsU3y3chWOez3Vjz3zfRr0Z2W8o5Hxm922HklGmY+nJvWLNRe4ycPg8bDxzDaZ2GK/H9a32R52zk2RrD3/8Y06cOQ0tP2ZAB5XuMwzSdFh9h0uDGcNYfzHLdGH3QYewPWHfwpLomghH403DUyZQF3FCr78uYOXcdgoJP43SINY+v/fUjjOpYUV11OXFFjb6v4MvFm7D/uMoXp44iaNWPmHB/bbVGR9R+8E1MnfYOHmtcsOuEiG4Pg4eHx79x7GwiIiKi28QRzSasweLR1RA97wWMOjYIs8e3R+nEMwiJNMG3qh88VAXNfP4vjO3/HDbUfA0/fvEMmrrGISzkEpzLV0FZN/0BLBs3AM/MC8sw44d93TWA4Fnof89U7LEPT+DYDBPWLMboatGY9/wIHHl4Nt7o7IPkC6E4E++BSlW84KLqleaIlXjtgRGYczrv4yq4NR6Fr398E519HYC0y4g4fQaxJi9UquwNV7WrKQdnomrPT+D9yHfYOOMeeFrC8N9HO+GtzTKOQhal7sc3Wz5HX08z9nzUA/0+P6aPL2/b+Ni6Dvux6mToj3um7kHGURpMle7FtP99ikG1VaXenIio0DBcSi+FCv5l4S5JG7UJU594GrP3Jtr+4kY8cd9XWzGrXylV9U6HOT4SIeEXkeLsg8p63wxITzyALx9/CB8EXlafd0STN1biz+dr624o2UnbNgXtH/oaYTc6BRnP5wujcGzQbIxvXxqJZ0IQafJFVT8PPV7F+b/Gov9zG1DztR/xxTNN4RoXhpBLzihfpSys2WgZxg14BvMybSznfOTYbALWLB6NatHz8PyII3h49hvo7JOMC6FnEO9RCVW8XFQ6mBGx8jU8MGIOcstGxlpj8dfqV9FYndbrpSN180S0GvQjIiwZrpv5Y/HyxVGYPbI+3CQwpFguzsOwlmOxSrKUqRLunfY/fDqoNlzVviRGhSLsUjpKVfBHWesJxqapT+Dp2XuR+Qx7ouW4n/DT2JbwNKptx51F6AXAp3J5eJoScOCnN7G08nS81tWALRNb4+Efztv+joiKmgzxWCIiIiK6dQxw7/Amvh7fGOd+fhY9WrRCp04BaPPwTOy8nA6TXx+M//gr/G/2M6h67FsM79wSbTp3RIu2g/H1vgTd4qLXuGfR1vo4O+8M7mj/1rd4o2U05o+/F62at0eXDs3Q9vGvsSderda3B8aP7QEP28dz5dkZE76QgIMJiQd/xuiuzdC8Uzd0bd8Ejds9ind+24/L7qXUBy2IWrUY/0Srmq6pPHr2b4PsnnmX7tYHHTyMSE/di6V/nbAGbfK8jTxwaoQXZs/AoNouSAr+HeN6t0Cz9l3RVaVB086j8MO+eBi82+O1T19CizylrSPc3Uy4tGcupo64B80aNEOHbj3Qrb3ax95vYGl4GuBaH0PHPoiKuqSdin3fjsYjAwfh8Y82QJ1qIG0vvhkxCIMGqdfAgXjo1d9wJpeAwFXqfHZ482uMb3wOPz/bAy1adUKngDZ4eOZOtW4T/PqMx8df/Q+zn6mKY98OR+eWbdC5Ywu0Hfw1rNmoF8Y929bWKiLvDO7t8da3b6Bl9HyMv7cVmrfvgg7N2uLxr/cgHib49hiPsT1yz0WW4z9j7CMDMWjIu1gZJQdtxok5L2OwLS3uf2UuIrKkRYkWY/HR8HowhazCzHEjMOz5Cfjsu9UI1sEXJzR6YTZmDKoNl6Rg/D6uN1o0a4+uXTugWdPOGPXDPsQbvNH+tU/xUpYT7NbxDXw2pgU8Dck4PncserYIQMeOAWjW6Wn8Z/sV1H1yBl7umHNbDyIqOhjYICIiIrotDCjhWxoR34/AY2//iWN6gEYLord9jmnzw1X1zoRKnXqh2tGZGPr4u1gRkqT/yhK1ETNm/IUIs8QHVIWtYT4rWoYS8CsTiR9HDMSrv+6DxBmkMhn593R8MD9Mb9e7Yw+0zNP8o0ZUfuQ5PFTFAbiyB7NenIglJ649A08M24jvxvVF6+7vWRfErMXCNRf0Nny790M7N+via7zQrXd7uBvTkbxzCZbox/353MYNGeHd7wUMb6w2nLAdn4x6Fb8fjL3a4iXx5FJMGvMl9iYDjjUfwRNdr9vBbFzE3OcC0KbfS5i17CCirjYNsSDu4C949/udSEs3wLVxazQtYXsnKhiBmzZi89EomCWwkR6HU9s3Y+PGjdi4aRN2nLi2T7lS59O3dAS+H/EY3v7zGKzZKBrbPp+G+eGSSSqhU69qODpzKB5/dwWs2ciCqI0zMOMvGaTThPLtuyL/2cgPZSJ/xIiBr+LXfdHW/TVH4u/pH2B+mGzXGx17tNSfvSHLRRzdtgkbNx/BhVSdGEgI321Li83YHZp15FIj3KtUhU/kMrz+2HB89NsyrFj0Iz75fClC1E4YvfvhheGN4aZSYvsno/Dq7wcRe+0EY+mkMfjSeoLxyBNd1edsjL4YMOIB+JuAtOP/xatvzcdR26aTQlbhwyeexPTtiTA52JqIEFGRxsAGERER0W1iDp2Dt6dtQkymWmwy9uzYh0T9JP8Qvn/jM+yUHgwZJAQF4YjUiI3l4O+fpwhEBmaE/joRUzdczFJ5TsLurbuRIKst448qZfJQLDSWR897WsDFkI4rO+Zh7tGMHT7sLEhMtAZl1J5j48KVOCv13rLd0b9j5sCB0bsr+rTzgDE9EduWLLO2Wsj3Nm6kJLr07YhSxnTEb/wVvx1TO5KF+fgyrDqslhtLo2lAfdvSG0uNi0Ni5sS0sSDy+ClEqzQ1uHjBR48hUdjU+ZzzNqZtisl8PpP3YMe+RFiz0fd447OdyJyNEhAUdEQHVozl/JH/bBSKXydOxYaLWQ48aTe27k5Q2zWijH8V28JCZj6HRW+/iQUh15+/kl36omMpI9LjN+LX346p1MnCfBzLVh1Wy40o3TQA9e0BHY/26BLgBoN659Di37Aray+khD2YNWY6tuoLk4iKOgY2iIiIiG6XK0m4kk2FOC02BvF6eTKSsqtIJcbicor61+AAV9f8D2J4JelKlqCGVdqlaOvTbYMr3PPSF8WxHhrUclCVQQvC9u1D1jpudpICF2KZVEiNZdC5Xxd42pZLMbRMt75o62FQldKtWLL8vHUfC7CNHDnVRZN6LnpdF6KuoHK9eqiX9VVH7ZFOcyP8Klaw/l2+GOFcqjwqV6+FOnXroo6fmzRCUGlqguMt6sWQ/flMQ2xMvHV5cpL1kLJIjL0MazZyhasMsJIvV5CUXeZV270UbW1xYnB1ty4qZOZTf+Cn1VHZHLMT6jappwdktVyIwpXK2ZzfenVUnrMGfIx+FVHBdk4ca9RDTflDcxQO7D19fUBEsUSeRURKNglJREUOAxtEREREd5o5XVe8cqbet9XqDIXYMj7dYrZVFo0w5mXFrmVtrRAsuBgZmW1l8DqpQVi4JBhp6UaU7tQf3exDYxi90KNPG7gbLIjb9CdWXbAdYEG2kROTH3y9ZV0mVBvyNZavXo3V172W4LW21tqug3NeR54womS9AXjl83lYvzsYxw7swNYN67F2zRos+7g/yuY0SugtZk7PpRKu3remskH9V1jUOs22c3er5ko1m3PIByb4+XrrCo2p2hB8vTy787saS15rC32GHZzhbBu01OjjAy9JhPQoREg/LyIq1m7R3YeIiIiI7jZGJyc42WrEaany7D8vzDi0aAn2p6XDWLID+vWwVkSN3j3Rp7UbDJYY/P3HGly0friA28iB0QGOJlmZBTHBG7Bi+XIsz+m1dB6++HmT9e9uyIhKA2bgryWz8NKDbVGzdArC9m3B+tUr9HpWBIYgiQ/5bxMjHBxNOkhjiQnGhhXZnNerr6WY98XP2BRv+0sHB+isofKGmXENomKPgQ0iIiIiyhNLUgISLVJrN8DN41qnktxYTizG4p3JSDd4oF2/nvA2GuHTszdauRlgubgOf66LtX2y4NvIliUBCVdkXRacX/khRg4fjuE5vUa+hOmrw61/dwNG7/swYcqDqFYiHdGBM/F4u6Zo1+dhDBk6TK9n9GcbbAO00q1nQULCFd3ayXJ+JT4cmc15vfoaiZemr0a47dyYExORJD8b3OF5S8ZCIaLbiVcxEREREeXNlXCc0V1GjKhUrVrepwy1hGPJ4q1ISDfArfW96FXODz17B8DVYEbkmj+wIeNEGAXdRnZSwnBaRi5V6yrnXwm2SUpuinNAN3QorYrQaUfx0+RPsP5MflqVsClH4UpB2OmzupuKDIhaKR8nOC08HOfldBjLorK/q3VhViYTTKwtERULvFSJiIiIKG9SDmDH3gRYVBGyTNvuaJltfdAIt1IlbT/bWXBh+SJsjLPA4NoK/QcPwb0BrjCYz2LV4s3INCFFgbeRDXMwAnfKOB1GuLfugQ55+JPcOHh46MEqkX4B589eH6hwK1MKzrauNNdJTYWe48XgBJf8jwFL1zEjOHAnIs0qR7i3Ro98nGBLyF4cuKD+0OCK5h3b4vqxc93Q8vnn0N2D1SWi4oBXKhERERHlUSzWLlihK5KmigPw6rPNVPUvA6MnGj3xBVasfde24BrLxdVYtO4SLAZnBDw7Gq1dDDCHrsDi7cm2T9gVfBvXS8bWX+chODUdprL34c2JfVAh24E9XVChQR345WHQz+SwMERIgxJTXbTv6JWhMG1E6YAx+HpKH+Q0c25aZASidGOU6qjfMPPUt1QwyVt/xbzgVKSbyuK+NyeiT/YnGC4VGqBOxhOcvA1L1khrDyO8+7yA0c0zng83NBr2JWY/1xBuOQWpiKhIYWCDiIiIiPIsbvUMTFsVoSqE7mg2Zg6W/e89vDR8KJ5+9g3MXLAOiz/oB/+Ua2NmXBOH9QtXI8JsgMnRCQ4GM04u/wNB2fTkKPg2rpe6dzYmfr0f8XBE9UFfYcXy7/DuyyPx5JAhGDpiDN6c9i0WbgrClqWfYXCt3CMbKTvnY+HhZKQbvdH7g9/xzdvP4amnnsH4GXOxcs44NIrchUNxEr24njl4EwKla4zRC/0mfoEJI4di+MvvY9YbvSC9W6gAUvdi9sSvsT8ecKw+CF+tWI7v3n0ZI58cgiFDR2DMm9Pw7cJNCNqyFJ8NroVrZzgRm7/8HOuiLDC4NMLzPy3El2+NxrDhY/D+Tyswb1JXuBzYhsMJ7D5EVBzwFkpERER0S1mQEBuDJLMZSTGxSFT12qzSE6MRk2KGOSEasVdsCzMyxyM6NhXmtDhEx+nODDY3WLclAbExSTCbkxATm4hsNgtLYgxiktV21fZj8zqVhyUMc8cMwetz9yM63Q01uj+FVya/j3ffeh4PtfRC9K5f8PqIj2wfzixx01wsPpmqR5pITwvG0j/2WbtmZJXfbVgSEWc/1rhElSoZXUbgtCF4fPIiHIxOR+n6vfH0uHfwwbRpeH/SeDw3pA8CKqXhxNol2BiaXSplkbIHn73wOuYdjkO6e230HvUm3ntvIl58qD7il72NQY9+hC0XVZomxSM+a2OU5EDMmjIHwQmAU+VuGP3O+5g87kn0buwLJ9tHcpTr+UxHYnQMUlReSIiORfbZSJ3nVDPS4qIRl5oxlXLOR5aEWMQkyfHEIDa7zKv+NjEmBsnqbxPVdvMsPQ6XolPU8STikjqu63Nf7teN3eXAaRjy+GQsOhiN9NL10fvpcXjng2mY9v4kjH9uCPoEVELaibVYsjE0U7pZQn/D2BEfYk1YMoylGuC+ZydgyuTxeLKLN0IXvIFHn/weh1MZ2CAqDgweHh68WomIiIgo31wrNEWHdo3h7+0CS9x5nAjaiE2HorIPVmheePDbDfisT0lY9s3APffOwJEbVFhF/rdxAy5+aNymLRpWL4cyLgYkX47C2dNHsDdoP0Jjc9mRrIylULNtZ7SuVx5uqZE4unUt/j4SnSWokj3nCs3RrXNTVCnjhJQLxxG4dj0OXCjQEVEmLvBr3AZtG1ZHuTIuMCRfRtTZ0ziyNwj7Q2OzDe5pzuXRtGtnNK9WGsb4cziyZS02HYuFxbk3/hP0DR4oZcaWia3x8A/nbX9AREUNAxtEREREdFsYfYfgfxs+RFc3M/Z83AP9Zh7NUyCA6I5gYIOo2GBXFCIiIiK6DYyo2G8A2roZkJ62DyuWHmdQg4iICgUDG0RERER06xlrYMD9LeBsSEfq3uVYepJhDSIiKhwMbBARERHRLWeqNwD9GjjAkJ6CoGVLcZpxDSoG2GefqHgwOTs7T7L9TERERER0SxiQiisXTmDb8h/w7dw9uMixMqmoS49BxLkInNy5CktWbkNoHMMcREUVBw8lIiIiIioETZs2xciRI7FixQr88ccftqVERHSrsSsKEREREVEhKFeuHPr37486derYlhAR0e3AwAYRERERERERFVsMbBARERERERFRscXABhEREREREREVWwxsEBEREREREVGxxcAGERERERERERVbDGwQERERERERUbHFwAYRERERERERFVsMbBARERERERFRscXABhEREREREREVWwxsEBEREREREVGxxcAGERERERERERVbDGwQERERERERUbHFwAYRERERERERFVsMbBARERERERFRscXABhEREREREREVWwxsEBEREREREVGxxcAGEREREVE+VKhQAX369LH9ljetW7dGw4YNbb8REVFhYmCDiIiIiCgfRo4ciW+//RY9evSwLbmxKlWq4LfffsOkSZNsS4iIqDAxsEFERERElA9//PGH/nfKlClwdnbWP9/I5MmT4ejoiEWLFtmWEBFRYWJgg4iIiIgoH4KCgrBgwQJUrlwZo0aNsi3NXteuXdG9e3fs378fv/76q20pEREVJgY2iIiIiIjy6f3330diYiJefPFFlC9f3rY0MycnJ7z77rv65wkTJsBiseifiYiocDGwQURERESUTxEREfj888/h4uKCiRMn2pZmNmLECD2+hrTu2Llzp20pEREVNgY2iIiIiIgK4JtvvkFoaCj69++PNm3a2JZa+fn5YezYsUhISMB7771nW0pERLcCAxtERERERAWQnJysBwYVErxwcHDQPwvpeuLq6opPP/0UkZGRtqVERHQrmJydnTnvFBERERFRAZw4cQIBAQFo1qwZvL29UbFiRb1cBgyV96TVBsfWICK6tQweHh7ptp+JiIiIiCifateujdWrV8NkMtmWWA0ePBh///237TciIrpV2BWFiIiIiOgmBAcH46effrL9ZrVy5UoGNYiIbhO22CAiIiIiukmlSpXCwYMHbb9BDyYqA4sSEdGtxxYbREREREQ3KSYmBjNnztQ/L126lEENIqLbiC02iIiIiIgKgYyx8eOPP2LUqFFISkqyLSUioluNgY1ioH9r2w88U0RERJRfBtu/dFtIcMNsNtt+o8L0Z6DtByKiLNgVhYiIiIiokDCoQUR0+7HFRjHAFhtERERUYHdTiw2TG9INTrZf6K6WngaDOUH9YLH+rrDFBhHlhIGNYsAe2DCq/0wmVToxGGAxW2BJT0e6eun3DOpdowEG9Uq3pMMir3TrF4FBfV5eJpNRfUmkw2y+9l5RI8dhcrA2JJJjNFuK5n7eLkajnHNrepjlnBeB9JD9kf3KSPJc2i16QmWUvOtg0j/L8Us6UO7s172kn1xH9nuFLBOyXO4X9t+FfEbS1/7Z28l+j8q4P4L3AaJCcBcFNswu1ZDu6KPvX3T30l9DafFwSDyiAxx2DGwQUU4Y2CgG7IGNGnWqoMeATnAu4Yxtfwdhd+ABXEm6ot+rVb86GgfUQ8Uq5RB++hz27wrGkX1H9XtSUajfrDa69++A5Csp+Gf5VhzaY32vqKmujrH/4Hv0z1vX7ETgP7v0z/9W9ZvWQs8Huuif//5rM3ZvO6B/vpN63t9Z56eMTgeHYdH/LbP9VrgqVPbDg0PvhcnRhD1bD2D9ss22d+hGSrg4o3q9qqhdvxrW/LEBcbHxernJaIRHKXfUaVQTtRpUQ2nvUtZAlaojnAuLxIIflyA19Voh8nZxcnbEQ0P7wbeij22J1c6Ne7FxFUuyRDflLgtslPariypVqtuW0N0oIvI8zpw+xMAGEeUZAxvFgD2wUa9xLVXB64uy5byxYWWgDlBEnovS77XrHoDu/TvCr4IPzoZHYOOKQPy9fIt+z93DDa27NceAwb0QcTYKf/22Gru37tfvFTUSnHnm9Sf1z0t+XYVl89fqnwuiVBlP1KxfVR9/yIkzOBkcYnun+GjTpTmeeOER/fPv3/5x9ZzeSY8/9xDadmuJ1JRU7Ni4B5djExBx5gK2rt9p+0Thql7HH2Mmj4SjowPWLd2Med//aXuHbqRcxbKo1bA6KviXw1+/r0Zc9GW9XAIILds1QaNW9VC1ZmW4ebrD3gDnZHAoZk76BinJqdYFhcjFtYQOpJRU1+XFyGgcDAq2vWMl57edylelvEvBq2xptGjfWC9ftegfLPr51gTNiP417rLAhm/lJqhRsx5279mPhETpqkB3EzeXErCYk5F0+RwDG0SUZwxsFAP2wEb1Wv66xUZtVVnZv/MINq3ZjmMHT+rm27L8nge7wNnZCYkJSfhnRSCW/rZKd0kpV6Es2nRrjk692+LwvuNY/9cmHFH/FkWFGdioXrcKHhp6L7x9vbB2yQasWLDe9k7xUZQDG/FxCfjsnW9xJuSc7Z1bo1LV8nh01ANwcDRhx8a9WL34b9s7dCNNAurDt2JZfT/YtWmv/le4urniiecfRs0GVdW9w4SLF6J1dw/pfnLm9Dn8+vVCpKYUfosN77Jl8Mjw/qhcrQIO7j6Cn79YYHvnenUa1sCYySP0zwxsEBWCuzCwUalydXz1zX8REXHB9g7dLXx9SqNCeR/4eTsxsEFEeZa5ozwVaSmqshEbE480VQlx9XDRzcmFR0kP9XLXQQ3h4uoCT083uLq7qLKMAW6erurzrrrvfExULK4kJevP3e3kCXElVYlyU8fuYBujgYqfc6HnMfuDH/GfKd9hw8qttqWUm7KqUCjdUcJOnsnUtUS6pZf2KqnvE2dCzuPHT+fgy/et6TvvhyVIS701Y6U4ODmgQpVy8CzlCZPJwbaUiIiIiOjmscVGMWBvsVHWzxuNW9VHt77tEX0xBts37NHdUarWroTOfdqiWdsmSE5KglMJZ+zfcQirFv+DcFWpadq2IVp1boEqNSpitVq2Z/sBnA+/9oTDs7QHylf0RdValVDKqxScnBx1RUiarh87dFJXLEVAp2bwKeeNS1Ex2KvWcTEi+mqFSQb99CjthobN68FfbSfmYiy2/ROEqIhLOuDi41tG9/f3reCjAw4yCOT5M1E4fTQUJ4NP65Yl8sQ4txYb0ozex9dLrauKWldZuNrWJV1sTsm6jljXFdCxCdp2C9BN34VU4ORptJDxBsJOn9U/34i0hPFS+12lRiXdnN+zlDsMRqN+8HXhXBRWLFyPpm0aokrNSnqUKzkf0jUoKdE67omo27gWqtasCKPJpFvYSLoIGfdEWpJI4KVyVbVudQ7kM7Lu2Og4rFjwNxITEnNtsSEVV0kP6XLjrfKHi/pdBoaVsRJOHglByImwq2nr6GiCm4eb+mw1lK/sh5Kl3dU2Ha4+yDt/JhIr1XZzG1g2txYbTVs3RBOVT4WMjSBpVlOd+zI+pZCWps57WIRuaRSr8lfl6hVUvqusux5IK4FTwaE4ERyCC+etXaxE2fLe6PVAF926QMaGkXxlV0sdi3TDElvX71LHEIHqdaqiUlUJaLnAYk5X24nV3R4izl3AlcQbB/Xc3F1xzwOdVdp44ozK94Hrd6JB0zooV8kXLm4lkJKSisNqHyR95RzWbVITvuV91DXnhMsqPQ6p7ZwPj0T85cxNo0uofCrXQNXalXVXMjc3V71c0iD05FmcPhaCmEtxOi9XVHmtQfM6+rjDVT49fjgEZ1Xeld/rNa2tzrMXItW5Cgrcr/JTnF5PRvYBeAc83hvmNDNWLvxb50lZt6R1j/4dUbtRDR34lG5EJw6fUufFjKCtB7B76z69DmldUaGKn86fEgSRdSYlJePShWgcO3BSn5/EDPm8ZBkPlKvohxp1q6jrxENfp/ZBQNf9tUkHVruoe5SM6+Ho5ICoyEs4qY5LBP4dhMN7M4/5U9AWG7Lt8upc1ahbVV9T0sonNcWs88AJdW+Q8xYXY+2SI2Q7bbq20D9v37gH0ereJuP8+FVQ15I632fV5+V+kRMZr8RZXXPV6vir67iCul+U1vnUfuzSpWeluk/I8QpZ7qPOn9wjpeWKdJOT6y1anUcJQAXvP46UK6mZrsEy3qVQUd0jJF+7q3S8HJeo8swZ7N95CC3bNUbFauXVveKK3o50D5N7TuvOzfTfbly1TeWfU/pnO//qFdH13vb656Ct+9W9/KD+Wcg9TvJ6zXrVdFrKPSPpSgrOqvvn0UMnEHXuks5HopX6PqjXpBZkWku5F0paVKvtr/a1gjr/Drii/i7yzAV97cWqNJd9syutvmckf1VR+bGk3o7kF+tAsSsWr0dEhu8nub9VVvvsX6MCPDzdoO5m6lq5rNLrrEqvE+qaVnm7iA6GXaTYb/R3AbbYuPuxxQYRFYTJ2dl5ku1nKqJqV7T+K4X0Ei4ldFcURydHHWA4duCErqhKYVwKojLWQQkXJ12gltYZ59XvtVThXcYpkADDzk179WdkEFE7+fuW7RurylRdXTGppAqR5SqVVRXvMno7CfFJ6pWIOo1romGLuvAp64XLMQmqcnPxamFVKlJValZGQKemqNOgOi6ril3I8TAd+PArXxbN2zVC44D6us9/1Zr+KO/vp8fAMJqMiFfrSkpM0i1KZIyQFu2b6HUeVYXWY4euFcqlAuGnKoXN1b5KAETS4eq6vOQpsElV1OJ1JU4K3fY++kIK7PK0WF77dhxC5LmLtneyJ8cthW8JXEgwqa469mq6wmxdh4Ojg6pg79b70FClmxy7jBsglVNp8i8VGHm1VMfSsEU9XdGTIFHC5UQ9loCHOlfSVUDWLZWDaur8SZcLWbebuxt2bNqNZFWRlGXyGSEVhNMqTe0kQCSVEDnORi1V2jaopgr/lXR6lFSVOqkpxEt6qPMnlREvnzKoryrMEnjIuk15SQBiu6pk5jYjhqS9VHglj0mlVNI8I9mfLn3b6XVKpUOCKLJNCW6Ur+yr86mDSl8PT3ddgW/Ush5qSLBF5Tk5xxJ8CM8QLJE88cCTfXVF8NLFWBzafW1shjqqgt774W56W/Fx8XBR14dsX45PKvEy8KgEkGRWDUn7WHV+bkRa99w3uJc+r+6qEiWVsxaq8iiVRQlg+an0lnPv6uai96u5ek+CBFJZ9PYro8+5VP4lYJGRVGbbdGmhj7e2uj78VX6RfZZrTAdg0iw4dyZSV+xKeZXUARDJ5xJgkWDKpQuxaNy6vtqXJmq7ZXEuPEJVNiMyBdHs5D4hA4JKJTPxcpIOBtnPqdwn+g7srvO3cC7hBL+KZfU5kuDgUXU/ERKYaqbyfr2mtVC9blW1r+X18UpwSmaokeCAHKMEPNw9XdU9pLr+fKMWdayVfJUeMpCxHONBdb4kqCOVafsMP5J+8p68Th8Ly5SvhaRL6y7N9c8njoTgyL5j+uec6P1Q507uL81aN9D3KdkPyacS6JA8IIEcaY0i+VX+lQpyXZVPej3URe+HnAMJ7Mqx11b3TLmmJdl2bNht28r15DhkHU0CGuhzK/dSCarZj12OI2jLPp1WV6/7Vg30q37T2iqP+uugqQT2JPgl16tcM3LOZUBXCbRJ3mvWtpE+JsnTZct5qc86IyEuUQcwGqjlHh5u2L4hSAeoGjSrrYNzsn0JIMpg0hlJAFvyuLwvAQt74EOCLNXU/b+52lajlnXVefdHRZV+fhV9dPAhJTlFXUNJSFTfBaJt95Zo3yNApa/KO+o7RdZnv7dUrq6uFfk7lX8luCZBcrkv6oFr1XUl91O5tzZoXlvfh+TaljQrp9Yl4z9duhCjW9lJS8RG6v7ZRJ1TOa4qKk9LXpQAkqSNXPMygHbG7zPKwV0U2Eh3LA33kup7rmQZ7Ny1BwkJ1jx5p0jZw8fHB7Vr10bjxur7uFEj1KpVC+XLl1ffSS5q/xL0d3Bu360Zubm5oWLFimjRogUaNmyI6tWro1SpUkhKSlLXedrVAOONuLq6omvXrmjSpAnq1auX7evixYt6/4oad3Vv9VT3JHdXE4yp8qDj2vEGh9t+ICLKgl1RihF5kidP/M2qAC5dUUqpQqNUpKTJecnSJXWB+PihUzoQId1T/FRFUb5wpYLmrgqBUui9EHFRP93LqIaqsEslSrq2pKrPXFFfnMJHFaBbtG+EmqrCLN1XpGJ9KSJaBxGkMC+Fa3sF3uTgoIMiUoGLU4VNaTEiT/RLlfFA7cY10K1/B/3UUbrGSCEkRRVEpUIlhVapuEuhPjelVGVdKpFSQZLKQ6Z1qcqerKupKgDLk1Yp6GZ8Op+qKuHxqiIgL2k1kBupJEkwqGOv1qjTsKZap4tuQRGvKseyjiRbGoaeOKMrmY7ODrrSIWktZN8k2CMVTHmaH3I8XFcMhOxf1RqV0a5ngK4MybakcG5ft1QA8lIAksBLvSa10f2+jjoQJRITk3R6SEW1UUBdXRlysnVRkiBIt34ddMVDAj2yP/ZtZjymwtREKruqkuTm6aIDIVKhlqfVnVS6druvA+qovCF5WfbZ0dlZV7ylFURByMxArTu3QI16VXW3h5SUFIntqOvDC606N0V1VSnKD2lZIa2jpMJmdDTp4IibqhxLxbdz37a6YiethmTaVFXUVuejtA48SKumrKRi36lPW30OZIpCCeSlXEnWFT+5llp0aAInR2uwIeLsBezctA8xl2L15+vbgioSNJRrUlqeSFAta0DJTs63VBQlL1y0tRSwM6sCsQQppSWL/t1suXr+7flTNFTnTGa+kfE45NpPVvlTzl1FVamUFjKSv4TRZEC5yr7quBupNG6mAwPyBF+2YV+vBBEkwJmo7kv2bC33Ivv7GbvJFJRca+XUPrVSadSqSwt13ZXUFWoJ/Mj+lPEpidZqecsOjXWa6mmzs2igjrmdOqcy6Kq6fPX1kJxN4Cgj95JuurWc5AkJbKWpY8l47BJMk2CV0Nd9LX+dfhKQkwCZnKNkWz6o36yWDjB7lS2jPy8BJNlXuTdLEECCbnKO5F5Tuaqch5bqnqoqTnm4d+aFtKaTwIRcQ+4l1XeBOha5L0mARYJbTVs1QLXalW2fvkaC0xLUkjwggVK5TiRvSyBdAn6derXV+VY4qH2XwITkleZtG+vvsLSU1AxplqDPl5AgvuT7tt1b6CCQ5D85n5IfPdX+yTIJOvr4XX+9Ed1OJUqU0MGMRx55BGPHjsXEiRPxxhtvYNSoUejdu7cOUMhn8kLKU3LNlStXDp07d8bzzz+v1zVu3Dg8+uijqFu3rirPuevP5UYCIW+99RamTp2Kd955J9tXtWrWVq1ERHcDttgoBuwtNtIlQq8qBvL0TirMMdFx2BO4H226tkQFfz9dEVo2b61+UikVKalIHgg6jMaqQCqVSSk0SnP9+NgE/bTSTrpbSEVr7Z8bsWX1dhzcfRSXo+NRpVYlVZh0Qvipcziw8wjiVeFTmhlLYVsq48EHTuiCu1Ti3T3ddVNzac0hTfjXL92sC6HSfFoK69J8fseGIKxbsgkbVwbixOHTenYE77Kl1R4YcGTvMV3xuVGLDWleHtCxqe5ysXPjbqz9azM2rNiK47Ku0h56XbIvh9W6pFuKTG8phV+xdd1OPX7A9n+CcFbtn1RAbkSe9Lfv0Qp+FX1xKTIaQVv34a9fV+luENv+3oV9Ow+pdJRgQJI+FzJQqTzxlObuUjF1UBVhaakiTySl4L9jwx7dokO2KwXxzn3a6MqhPMmUpuBLf1uj91HWvSfwAGIvxsKiaoE3arFRp0kNnRe8fEpjkzpv61V6bFmzXXfHkebb9krmwT1HdaCnUrXyOv3kybGk0d/LtmL1H/+oNNmtj2v/rsP6fOYmtxYbtRpU111ERNT5S9iz7aDOWzJgrZd3af3UX78XcQnb1gepfd+hZ+OQ1jCShtIdYMemPfozQpriy9N7edp/6lhYphYbUoFvpPZHSEFPBsKUcTgkX8iUx/JEV4IlTiofh54Ix9GDJ/VncyIBttadm1tbayQm67y8dslGBG3eq1uLSKsTCRbKU2vpNrN8/jocUNuRdK6s0lcCKgdVOmZ9Si6tcKRLx6ZVgWrftmHX5j0IU9eVBE/kOpBgozxxlxZQ6WbpOmTRrVdKlympBwCV1gNyTUve3b/9kH4SL0GJ7MjT+6ZtGqi8Fa/3Q1oR2UlLK7n2JGgnFeUQlZ4/fPqrvi4k70orISFBAun+Itex7s5w6LSusMo5cnB01C0oQlR6yv2hXbcA3eJGXXp6v2Rw13+WbVF5Y5de75nT5xF5PgrhIed11w9pvXRI3WN+/mKefl/2QfJSRvltsSFdstr3DNDnWtLl8J5jerBgyQtyvBJE8Cjppj7nop90SncnCa5IxVuCfyIuJg5HD5xU18VmbFq5DYHqXinXnFS6c1KqdEm0vydA5Vt3fW7+WR6oByrW17I6tp0q38ixyz7Jue7Uu7VOW8k7a/7coLvpyH1UqigSrJW0PXn0tG5R5ubqooNnckzSnezYwRNYsfBvBKp1nwuN0C0c5N4tLfSkJZKktxyTBAMkYCr2qmsva16UgLK9NZt05bC32JBub9JVSPLyWrVvcn8I2rwfkWejdDcQSWPJS/I3QgJfEszV9Sv1PzmPG1YEqjyzSe3rKX3fk9Yozmr/ZLBr6aIlea5Dj9b6u8Ws8rrch6TbnZwnufdJdz5pZShBqTLqfi5BETlO6aaz+g+VXupaDFb3EX29qXMn51WCxrJuykXu9eBio6i12ChZsiQGDBiA8PBw/Pnnn/j111+xePFihISE6JYcPXr0wJEjR3DpUuZAc3YkqCEtLR544AHd6uOPP/7A77//jq1breNLPfjggzh+/LhuaSEtN27Ew8MDDz30EObPn6+DGwsXLrzudfToUf0QoKhhiw0iKgi22ChGpHArFQCpsEsFyEMVEitVqaArtvKE/+KFGN2XWwYYlYKmFAyl8ixN/1NTU3VFUlouZO2PLNOgSgVTKiVSWTlxOASnj4Xqp7pGo0E/kZa/kRYh4SFndRcUKaDK2BPypNHV0xXV6lTWFfy4WOn7fEa3OpDCvAQcpHmxrMevgq9uLSBPH+s1qalbYDiXcFYV15L6iWtupM+8rMuk1iVPymUdurm8rEtVDqUiKJVmacIcrSqhsp92Mm6F9EuXV3bN97OS45PxH+Tpt6TFri37dFBBCtGyjogzUTpNZF0y3oZUNKQyLJVwaWotT1vLVvBBCXV8UvGXVh32bjsurs668iHjYcgYCjs27VWVmVCEHrfun3w2zfbU8kakxYa0VBHS/aRes1q6ybY085YnybIdCW7Ym//LE1E5v9LMvYzKMxKMkc9K6w0ZH0IqMPlpKpsX0g1CKleSx44eOK7zh5CBcPds2a8qjkf0GA8njpzST66lBYDJsWADvUpwZOemPTgUdFQHtk4cOqX74QsJNMk5yQ8JaEkFVSqeEmALP2Vdl1TkzqrzLcEXqXTLdSPdOISktRRMs7p44ZIeb0QqyqdUfpLuFxKgkbE4ZH1y3HaSr+T62b/jsA6syLUjrUDkepS/C1bbkwpsdudKgjuOJRx1AFIChdItJiOpAEhFN9nWOiMp6crV6yJjN52Du46oSvk+HdAIO3FG/XtCH7OQ1g6yz0I/Wazkq69zqfRuWbtDpddJ3ZLJvl5plSBjeZwLi9AVViEBNPv7kvdullzzEoz0LO2O6KhYbP17p6qwW/dDAgJb1u7U9wTJ67rFhvH6vCBBjW0bgnBEVdzlWpd9k6DEjUiFIERdtwnqPOquHKrCLi1wpOuLvCdpbQ/aSIsNuWfKdS8tSqSLoLSQqN+0jh6nRfKn3AtljCQh14G0yJCWdHIeJYgo3dnkXiF5UvKfBO4Ki9xTJfgm3ZMkYNGoRV3d2kpaacgYGBK8ke5jWUm+lHwt43XogJe6R8o1IeMNST6RMVUkyCEcHRx0yz25R15Q31Vb/96FY7bzJOkddurM1eCaBFJkXChpLSfds6R7mNzv6zerowPO0q1F9lc+R3QnSfeQzZs3Y8uWLdi9ezf279+PQ4cO6Z/PnTuHgIAA3coiL5ydnXV3FunCcvr0aWzfvh379u3T6zpw4ACioqJ0txRpBZJXZ8+e1fuU3Ss+vvDuIUREdxoDG8WIVGSkWbN0R5GnyTIQpIyrIBV+aTJ8SRUUpVmv9H+XJ/TSP7lG/aq6MC/NqqWgKRUiO6kESYVAAhbyeQmCyNNBefIs68yuqWO4qihKIVQqNNKCQirI0pxYWkZIf3N52iYFbzsJNti7Z7i6u+rmzpVVYdXLz0s/AZcK7VlV4UlNzr05ugRB9LrUbkngwb4uaQ0ilRapIJ8NO4/CmKpSxoCwzt5g0C0wQlVhPbvKpLSqiIqIVpWoEF149ynvpdPEwWTSrU+k24EM8CiVPnvaO6jCvXQdkif/Emw6pSr9EvTI7wB4UqCXdBBSqZInsZWqV0AZ3zK4oNYrT+DPh0tl0poe0i0hSFVWJVAjzdcbq0pLt/4d0apzc9RW+UgCZHlp3pofMjCs5BdpPi6VW3tLGQliHNwTrPvmS3BIWhHZm6AXlByvtH6R4J7kAb3euII/yYu+FIt9an1S4ZdAg4xVYid5Qp4cy3WXEJ+gr6GcSAVMuo/JuAGe6lqpWrOSqgD7664LMlBvdmT/Jd0kf0jFWM6KBCHOhUfqc5oTeYIt3WVcVSVa9lnSPD8kiCkDRsr5EPK03D7YrLTSykr2Syqpcm1K5V5aB8XHxOvWHbeV2m/pPuGkKgWXVWVfzo2MByGtM2SwTRmcNCEuQQcq5VrJ7un1qaNh+u/k/plT4CgrCR7JmBAS+DJbzLrFQ6debXQ3ERl0VgZ8lRZS4up1r+650s1CWqXIeBe+FbyRnJKigxahp87qcVGE3GMlcCxpK/ukj0GdF7lXSF6S7UqQOz8kLxqzCeoIOb+yX/KdIC1R5NxLNywZ+FUCPXLfuHTh+ifOFnXccp0c2Xtc3+ck3SRgl10AWY5JupHIeZBjkgCNpGF24wXI/VECZrI/EmyRcUhkbBrfij7q/KThqKTXSWtQiehOkjEq1q9frwMPMTEx+hqQlwQhIiMjdVBDrv+8kC4rMq6GBBx27dqF8+fP6yBpXFwcTpw4gU2bNsHf3x9VqlSx/QUREdkxsFHMSBeSy6riIIEIKXDK+BdSSZVuEVGR1qezUhGTwr00UZZB8KQ1g/TNjr102dqdxUa6bUgFqHFAAzwyrD+en/j01dfD6nd5yp2VzIAhT3ClEF+jnnVMDXmaJoV4KayeCzuvKgjW2Q7s7NWDNUs24utpP2HWlO+vvd79Af/9bI4eCDUv7F1o1v6V/bp+mDHn6gwEN8u+3xKckFdOpPIpzbmlsC3dT6TJudHBqFuoyEwSkeeuH7H96rpVhcz+NLMg7HWvFQvWZpsev8xeoCt2QvZTnvoun7dWz2oTp/KRVGaatmmkxy2RcR4kTxQmidVk7PaUkex79u8UjFT4snZpuFkSuMrJjd7LSCqmVWpWxH1DemHk+MevXmMvTBymA1E5kZYa0gLEYMsfEvhJN8sAdPrXbEnrLOmeIMHCvHQrykq25+bppsc1GDL6wQz3hGHo3r+j7VPZk5xjbx1U5OTxXOWXBI/2bN2vZ6AKXLdLj2Ek3Z5koGUZU6Vr3w76nNjZgyXS3ezbj3/Gf9Q1mvGa/eK9H64O4CokD0isUVrrpam8nfEw8jJ4YFaOzo5wVq/syLpl9XKf+GnWPHz5wf8y7dsstW9ynNmRMXJyaxafVV7yixy3OBB0RKXX/+l9yLhPkl4Hgq51TSMqSqTVhQQ0pItKcnLeptmX1lEyqOeVK1cQGnrtIZG4fPmybr3h5+enX0RElBkDG8WMVG4kgCFPiSWgITOYSKsICWRIgVQKgjKWg4zdIP21ZUYAaUkhLTbk/YwtNmQQuo73tIZMjSpPemV8jv/7Yj4W/7ICMvinfYBBaRItQQt5mit/H3Xhkp6NRSr70lRZZv5wdnVGyEnp63wBqVeuVS7lKam0JJFCc/0mtfSgedLSQFqeyOB0+pXXp6OX1brUMUhlWWYKkKejeV2Xq22Mg7y2SJDm0FHnonTT+Rr1qulBSaWvuDzxzCpFFVhkXAhpzSL9+KUliRROZNDKK6riI91WMpKnL5FnLyAlNVXPYiIDTso5ym+lUGZQiI6yBnFk5hWZolEq0dI6wJ4eGZ88y7/SFSn89HlsXLkNC/+3FH/83wrEXIrRTedlPI88Jk+eyTazPbdqUY7vFZAEUApzfWplOe6jLMrr/tdvXhuderfVYygEbdmP37/7E3O+XqTHFrFP/ytPsmWQWXsekLwm4yTIU3PrFL7p+km3zJaS3QCOdtJ6yrecD86cjtD3iPymR2mf0ujzUHd1XdfXAdQ/56zCj5/9huXz1+ouA3aSv6WLgUWl+cUouR8l6pZKXXq3g4+0Usjl6aRusVC2jPpcwbodZSUBW7nvyQC/0vJIpiOWQY7lab902WrbPUB3y5JuMTJVdrqtwpxRXs9nRvJ5aZ0i17N0DVn4v2VY8N+/EHI0TLc2kJlF5FwKCUbLdMQSfJNubs3aNLraCubq/Uu97C235F/pxiafl7Rq000dg3dJuLq66O4cvR7sqtP8RjzLeOr7ntxfZHDnrn3bqfVYp7fNKlqdx/i4y7q7SbvuLVDB3zoDUMZ9y25cF0myvKadfH9JazIZ00kGvW7fvZU6Bmv3wawkUBlpm8Grgn95HXiV2VkkIpJxnwoS4CG6VawtYR30oJwytkaNGjXw3Xff6e4geSH31goVKuhrJTo6c1dCCY5ERESgdOnSemDQvHrssccwe/ZsfPzxx3juuefQrl07lCljHaSYiOhuwsBGMSOFx9hoa1cUqQRJwVj+lUqIjCkhZUupLMnv8rRPun9I02dpFixdETIGNqSfeZVa/pApOKVCLl1XZCBDafYrBVizOU1t0Dp+gwwAKhUv2X70hWgc3nscaSlpqFytoh6IULZ18kio7u4iBX17IffC+Uu6i4gECKR/uRROO9zTCh17tUFn26ttlxZ66r7cyLqkZYR0rZB1tdTrao2Ova+tq52q0Oim5ooUeKVQLPUEqUi079lKNxOXik9upEIhzculxYX03W/ZoSk63tPm6rakb7w9yCFpJU2vpV+5Ti+VntKdR/rGx0RfxoUsU8vKrDXSF12anMsYGTLzTAe9b231uuW45BzkRrr9nFYVKElqmdpTAlSStjJAoT09ZLBQCXwJCbh0VO/L7AUyCKZULOScS7cZCWjcqFUKFZyMXVClVmWd1tItTCq80ipCWj3pVhgqf0qebasqnDKWhlTyPEt6oHHLerZuVjHYv/MQZGwU6XolYw1I/sguyCZdF6QV1ZnTZ/V0oPkl40BIKy9pdeSu9lWCKVIxNqkKrkyHag92Vlf5p36T2up3M0KOhesxIOT4mrZtqK/BjNe3DGQrVNVXVdJT9Dqki4YMziufkwEmb5YENGUQYRkXRdJXBsdsZ1u/DCoq9x3pXifvy/grkvaFQbrGyfG2UtdZvUY1dSBZ7qHSMkIuKplC2d4KSoLPB/WAtgn62mum0sp+T7Knlbyka4owp5oRduqs7rJRyttTz+gix9KhV2t9j5ABZWUcnayksm+fgUYG8ZVtyH2ybdcW6rxV013ssiPjt5wPi9QD6NZvVlcP3irbyrhvElC+GRKskMGNpfWgBJokuCOzxNjziwTaJf2E3B8P7A7W4/J4+5XW42vIZ7Oml6Q5UVHQtm1bPPnkk/olg39KgCIsLEx3H5GuJHkhQW4Z9FNkbeUhZZrExETdXUXG4siNdJFZtGiRHqdDusRI9xZZt3R16dKli555JS/rISIqLhjYKGbkSaNUumWMAjuzqihIX2vpdy2VBykIS0U7Ixnj4FLUJV0RyUgCHTJQZckyJdGlbzvcO6innoZPAg0y4KYURGWazO73ddIFXnkaEasq6zKtrPzr5uGi/1YCADIo4sWIzE8YZDwAmZ0i6vxF/eRSBn7r9+g9eGR4fwwcOUC/7hvSW1c6chNxRq1rZ7AOcEjFrl7T2mpdPXU3Gvu6Bjze+2rBWFpdyBgOMn2tDIbX5+FuGDjiPtinRr2RGKlM7jps3W+1rTqqQP/AE32ubqv7fR1gzPCUUbZ1eN9xPX2qjHUhMy0YDUb9FDnrAI7xMQm6T7p99pSaDarrbgoPD+un193n4e55mv5WugXpKXjV+uWJbENVEe7/mKTtfVfTo+8j3fUTXlG5agU88GRfPGx7/4En+6DLve11ZVaCYjK4oy0eRYUoXV1fEtiTtJXKWa+HuqJr3/a6lUxycqp+el2yjAfuHdhDV/Kle5jMnFFbVZTl6Z0MaimzysjMD17epfSYFxJ4kBYTGUlAxL2Uqw6enD8bqWfsyS+5v0hA06wqxxUql9NjsPQf3At1G9VWlXUZFyFOB/JkelOZstOijuuY2j/pniZdEWQsDvmbh56+92oerFDF2mRaKuoR56J06w4ZLPOeBzurz/XTQbabJWMLHQo6ogeulFZrUunvrdL5YbX+Xg90RVX1u7Tsksq7zMAkx1cYJIjSZ2B3PCT3hRED9DXc68EuOogoMxuFHg/T51jIfVkGIT4bcl7/br8nPfjUtbSSlwRthQSyZAYSGbRVKjvV6lTVxyL3vCZtGugWWxJsykoCJ/I38h0hgYg+D3XTfyODNsv9P6cxWiSfHTt8GonqHiaDmMpsW/ere3PGfZOZU26GBJoP7zmKM6fP6eB35RoV1DF10edJ1v/IsPt0YE5cvhyvWzjJAMsS/JNuWz0f6IwHh2ZOLwmyERUFHTt2xMiRI/H000/jnnvugQzavmfPHj3uRn66akmLDylrZdcKSpbJe3lpfRobG4tPP/0UEyZMwPvvv48vvvgCgYGB8Pb2Rs+eqqzXvDk8PbMPdBIRFUcMbBQz8rRTZq/IOCjgZVXZkGCHBCGEtNiQQIetRbOuiEjFOyFOBmm79kWZmpqipw6UqRHl6a6sW8bIkCkYZRrLDSu3IVK6fmT5bpVgiDx9lBYKEkyRJwhSsbmoKscywGFG0s//+KGTWPC/pXp2DKmESyCmIGQwwBNH1Lp++gt7th3IdV3S5Hn14n9w6li4foKZH3pbqpC/5o+NOsBxOfrGT1tk8NMwVQCXQIaLq4uufErXn8uxl68bi0Gaw8s+/b10iz4OGRMl43nJK6nAHNl/HHO//0Pvowx2eaP1yOwTJ4+F6nMtFW35bGpKih40ceu6XXpmHCp8h1RFetvfQYiLjtUFXWnxdECdr6VzVmH9X5v01KNX87Eqq5arXA6tOjaBq7uLbtmzc9M+XRndE3hQ5alolC3no5/ASzeDjEqVKQknB0c9WG1yUirMBchTMmX06kVyzYTq61ruG6Enw7Fp1VasWrgO2/7Zra9pe5aWTZwNP6/Hbtm4cru6B8jArdmPcyJ/98/yLXqATrkuCpPc4yQ4t3ntDqxbskHfx2TfJV1lxgKZilWWb12/UwcUC5I22Um5kqpbTcngxamq4iLpIgPknlPbk9mO1izeqINB+rPqnEhwauWi9Xr2mHNnLmRqQZeV3FOC1fUt9+i96tzLPShN5R9J44M7g/VUzTITVlYhx8P0NL0S5JGuOalpqSrvxWHP9kN6qu296t/sSMsgmVXor3lrrPeJDOe5sEhrkrDT57B5zQ51v9mpW9BIumVHjl/GJ5J78MaVW3E+LEKlV+GOoUNUmKTLyeDBg/HUU0/hgw8+0AGK++67T89iIuNt5JV0V7UHMLKSZdJyI7ugx43o757oaOzYsUNPHysBjiFDhnCsDiK6qxg8PDwKuehCha1/a9sP6kzJl5qMdVGtjr/uIiGktYRMRSlTQwr5TAV/meHEXzeDloq1FKiPSleJDKRFgXSX8KvgjTI+ZfSYG5djEhBxJlJVcC7Ds7QbyqtKloenqyogX8G+nYf1IHHyhSoj/cuI+T4ys4OjE2JU4f1UcJguSGclXWWksu9X0UdPCSstBKSBtp30Id+745B+kiszBTRoWlsvP60K5jIaf0b2dcnI+DK1atZ1ScFZCu7SYkX20U3tuzwFL1nSQzcPl0/u33VEVQgyt6LIjh7Bv7QHyvp56f2W8TzsW5IpdaXVhX32B0lLSb8HnrxXP5WXD+7esl9VpHbowEFWco5kUEE5Xm/f0rq7iP04EhKSsE+lhzzdlDSr06CGXi5PU6WVRkbylF7GWpEpLGUQ1xKu19YjEhOv6LSVSp6MOVCpSnm4e7rDwUmeCFkDZdIyRyp70aqSlHFMjpw8/txDegwDCax89s631+2TTGVbpbp1KrqsaS3T/Up6Sgsi2S/7zAmOTo6664V0rbqkKlj7dhzWy4U8FW8UUB8mtcPhIed1dyS7chXLonaD6vrno4dO6QpsRpJHZNpKEXLyjO6GcCMyRknjgLq6lYsMQptxYELpkiH7KGT2CpnC1k6uRXvLgyMHTujrzU5aP8mMDnKOnFQe1K2nVFpLNwPp6iFjUsh1JOdNgiASdJJpPiV/nFfXYoRal1wb0k1Fjlfyv3RHO30sXC+3q9e0Fiqq8ysBK5laVs55diStG6njkOmis6a1XDOyT+Uq+6KkdFkwGnTQVFoZyFgy0n1A7i0y5o6kz6HdR/XfSYshSWu5LmX8DLlO7bnwoPqMDKopy2QWn3IVfXVaSp6Xz2SXr+s0rIExk0fon1ct+geLfl6mf86NXAsS4JFWUxIYkmmkpaWIBDal65YEbjLO1iHdGKS7hpBpamWq5fyQ661K9Ur6PiF5R3oHyS1BWlNI+sg2Zfv2cTPkuneR8UXU+ZYWO3J/kvv5tSsWOLz/WKZxeeS6lTE2vGRMEsk/ulthtG71IIO61qxfVXcH+mzSNzpv6fueSmcZ00W6GEpXROkCJHlOrnNpel5f5RVx6liYngLaTvKndIOSVm1yn3BS97SMexdx7qKenUXIrFgVVV6QING+HQd1C76MKlergGq1rOPBSL6Wbop2MpOOTEcuLY90frF1g5PvKrkvZJx6WM6pXD8+ft46vaQLVsb0ku81GSyXcpEx0Yo5s0s1+FZugkqVq+Orb/6LiIjrB+i+U+Qal5eMgSHjbNStWxe1a9fG999/j1Onrn135US6h3z22WdYsmSJfkmrCzu5dmV8jcmTJ2Pnzp349ttvbe/kj5ubG9q3b48PP/wQzzzzDLZt22Z7p+jw9SmNCuV94OftBIfEI6r8ey0A+meg7QcioiwY2CgGMgY2qOiRQowEPWQg0K73tkMtVdGWAMuSOStxKOhooc3SUlRIYEOaqUsXiiVzVulKllR+pZk/3Tnd7m0PJ1W5lsrqiUMndR4sTiTwUbWWvw5WyqCp0mVNSKurvAY2/k3kvjNm0gg95krGwAZRthjYuO1kLIwOHTpgypQpePbZZ7F7927bOzmTbiKvvvoqgoOD9dgcR49aA4lCghoSKOnbt6+eCnbhwoW2d/JHZm5q06YNpk+frgcTlTE4ihoGNoioINgVhegmyRNX6RbQtV97VK3tr8fkkIFVw0+f012C7japKWbdf19atNw3uBeefulR9Ly/s+1dulNkmleVFXE25CxkZpziRlqSyFgQ9vwkeUxeOXVVICIqCiTIKN1O5DtRXvZWGzKuhgwQLeNYyPt5Id1QDh06pLuuVKlSJdP6ZODPJk2a6JlRzp/P3DrRZDLpMZnkc3bysyzPuA55SWsSV1dXHDhwQA8wSkR0tzA5OztPsv1MRVRta6t+KqJkmtSADk1Rp3EN/bRZnppKH/9dm/chPjZBD+h6N4k8H4X9Ow9jx8Y9V18yiKm02qA7R7p7hJ06Zx1UMp/9r4sCGbhUgoG7Aw9kylvHD55EYgEGQr3bSQVFZj2SKV9ljA8ZG+lGY3bQv9xd1GIj3bE03Ev6oWTJMti5a4+qnOd/BqjCJIGCoUOH6hYVEiiQl3QbadasmX7JbCT//PMPLl7MPENadqQrqMygIgGMqlWr4syZM3p9Mj2rDPY5YMAArFu3Dvv379fjINn17t0bo0eP1t1Krlyxdrdzd3fHww8/rIMYsk4Jmvj6+qJbt26oX78+1qxZo1uEZFxPUeHu5gJPDze4u5pgTJVubNeC9cHXes8REWXCrijFALuiFG0yin+VmhX1TA9Gk7U/e3RUDPbvOKIHCiUiKmwS2GjTtYUep0cGjN66bmex635EtxG7otwyEtgYOHCgblEhg3RKKw15SUBBpmyV4IGMiXHpUu7dUqV1hXQVkVlLJPggfy/rtC+XVhnz58/X43XIcrthw4bpwUBlP2RqVyH7I79LwEU+K4ENGXhUWoPITC3Lli3TA4rK8qKGXVGIqCDYYqMYYIuNok0GJpTpFWUQO5me8bgeDPE8KxlEdEvJALRyz5EBQKXCQpQjtti4pSQAUblyZd1CQ1pbVKhQQbeSkOle165dq2dnygtpsSFdWCTgIN1IZCyMhg0b6rE3pPuJBDXOnj2rt5eRv78/vLy8sHr16qvbkmCItByRcTlatmyJWrVq6XE/ZPwOmRUlNDRUB2CKIrbYIKKCYIuNYoAtNoiIiKjA2GLjlpHWUzKGhrSokJcEFCRAYW8lYW8RkduMYxnZW2hIIEICHBK4lPVJ0ELWkzWQKa0wJIghwRT7e7Jf0sJDltvH35BAhn2fJIBSVLHFBhEVBAcPJSIiIiIqAHsQQ8bCkJYWMpaGdDu5fPmyblkh7+cnqCEkOCFjZUiXEVmfrFfG6pDARHatsyTgIZ/N+J5sUwIYsh+yP7Ie+YyMqVGUgxpERAXFwAYRERERERERFVsMbBARERFR8WKQ6UytU6zydXe9pNsMEVF+cYyNYoBjbBAREVGB3YVjbPhXqYn58xflaRpVKl4cHE3wLlOKY2wQUb4wsFEMMLBBREREBXaXBTZcS1fVs6IcPXoMV65wWvW7kUuJEnArYWZgg4jyjIGNYoCBDSIiIiqwuyiwYXEuD4vJ0/Yb9bp3ENavXnTdFLB3A0N6MkxXQlX599q0tAxsEFFOGNgoBhjYICIiogK7q4Ys4PBwdlVrNcX4qX8icP08/G/Wy7aldxsp/F4rADOwQUQ54bcDERERERUTMqUpX/IaO+k39S/QusvD8C5b4eryu+vFp3pElDdssVEMsMUGEQmZAcC9pDvKVfRF6IkzSExItL1Dd4uK/uXUN7MBEWcikJZmQXo6b/xUCO6iFhvpDp5IN7rYfvv3at6qO0aO/cj2G7Bnx3rMnjHO9tvdwZCeBkNqtPpJAhxWbLFBRDlhYKMYYGCDiISTsyMqV6+I9j1aYfn8tarye8H2Dt0tuvRpB6PJiC1rdyA5KQWW9GsFeqICu8sGD7U4eNt++3dydHTCe1Nnwcvbx7bEasbH7+DQwX2234o/gyWRg4cSUZ4xsEFEREREVEyMGTMG48ePR3h4OCpWrIizZ8+ifPnyOHbsGLp37460tGuBACKifwuOsUFEREREVAz4+fnh+eefx5UrVzB58mS9LDAwEIsXL0bNmjXx9NNP62VERP82DGwQERERERUDEyZMgKurK2bNmoUzZ87YlgLvvvsuEhMTMW7cOPj4ZO6iQkT0b8DABhERERFREdeiRQvcf//9CAsLw5dffmlbanX+/Hl8/vnncHd3x1tvvWVbSkT078HABhERERFREWYwGHSrDCFdUJKTk/XPGX311Vc4ffo0HnroITRr1sy2lIjo34GBDSIiIiKiImzgwIFo1KgRNmzYgOXLl9uWZpaamoq3335bB0Hee+89GI0s5hPRvwfveERERERERZR0L3njjTf0bCcTJ060Lc3e2rVr9atx48YYNGiQbSkR0d2PgQ0iIiIioiLq4Ycfhre3N77//nscP37ctjRn0mojJSUFzz33nG0JEdHdj4ENIiIiIqIiSrqf/PDDD/jkk09sS25MxtmQ2VG+++472xIiorufwcPDI932MxERERERFQPS3WTZsmVYuHAhXnjhBdtSIqJ/J7bYICIiIiIiIqJii4ENIiIiIiIiIiq2GNggIiIiIiIiomKLgQ0iIiIiIiIiKrYY2CAiIiIiIiKiYouBDSIiIiIiIiIqthjYICIiIiIiIqJii4ENIiIiIiIiIiq2GNggIiIiIiIiomKLgQ0iIiIiIiIiKrYY2CAiIiIiIiKiYouBDSIiIiIiIiIqthjYICIiIiIiIqJii4ENIiIiIiIiIiq2GNggIiIiIiIiomKLgQ0iIiIiIiIiKrZMzs7Ok2w/ExERERFRMWA0GpGeno7AwEAcOXLEtpSI6N/J4OHhkW77mYiIiIiIiIioWGFXFCIiIiIiIiIqthjYICIiIiIiIqJii4ENIiIiKqI8UL//KIwZMxL31nOzLaNiya0O7n9jNhasDUTQnl3Yum4hZg1vanuT7n5uqHP/G5i9YC0Cg/Zg19Z1WDhrOJq62t4mGyP8B3+DjXv2YN27XWzLihCP+ug/agzGPNMfDTxsy6jw2dN5VH/UL2rpbKqPF+bvwt5t8/FSMyfbwqKBgQ0iIiIqmhxr4v5xb2H8+Al4+b6acLQtpmLG6I8hs+di5nP90bp2OXg4OqFU1Zbo90BH2weKIyO8Wz+J16dMwetPtII3S9Q3oCrrQ2Zj7szn0L91bZTzcIRTqapo2e8BdKzIhMvEJQDDRvdCNS8DDm/fbVt4A04l4VvBF5636eboWPN+jHtrPMa/ORb9qudzo0ZvtH7ydUyZ8jqeaOV9c5VQoxt8KlVCxbKeefhecERJ3/LwcTPZfi/6rqbzW2PRP7/pfKuZjyFwXwJKV2yFp0f3K1L3Pt5NiIiIiG6GsRT8GzVDk5reKFrPr4oGx5ZP4ZnOXjBZIrB8fFc0qN8Q9eu3Qf/xC2yfKIaMFfDAq5PxwrBheH7KOPTzYZE6R44t8dQzneFlsiBi+Xh0bVAfDevXR5v+47HglMX2IRKlez6JAf4mmMP/wq+rYmxLs+OKGv3ewYJt+7Br+05s+aQXnG3vFNgtvo8ZKzyAVye/gGHDnseUcf1Q4EvGWAuj5wchaGsgtgXtxZbvH0fNG8QsSg/8Djt37cCufz6wLaGbk4Jdc+ZjX6oBpbsOwcNVis69j3dhIiIiopvg3GMylixdgr9WzcQjXraFZGNE5ZbNUUlVPMxR6/DL/BNIVkstieHYeyDc+pHiyBKPyKjLsCAdlpjzOB/PCnpOjJVbork1A2DdL/NxwpoBEL73AMJTrZ8hxeiHex/phjJGM04tW4jtkk7ZMJZuhqGfq/vNlyPQuqwDDOo/k6PDTVfqbvV9zBIfiajL6jpJtyDm/HkU+JJx8IV/JRcYDepngxPK93wTUx73z+H4jTA5O8NRfdbgeNOhH7KxnFiMP3alqEzTDA8+1ABFpS0MAxtEREREN8OoCs/qH4PBpArbLFplZoJPWS9V9VL1majziDBblxZ/0Vj8fAe07dwFbTu9hOUJtsV0HZNPWXhZMwDO3z0ZoNAZK/RF/9auMJhDsG7lPlwf83FExa6v4KcV8/Heg7XhlhiB8xIoKCy3+j4WvRjPd2iLzl3aotNLy3Hzl4wF8fHxsBg90e6FcejlxXvvbWMJx+o1+5AGB9Tq0x+Ni0hvGeYAIiKim2RUBcIbMqqCou3HPFPrdDTd4a9ptQ+3fA8KdRuqYJ7buciPO3wO7o58ZYCTs61huzlNFYTzSuWLPOymUR1LQY4m17TNi9QYhB07hvDY3CrreTuW/LgVecNodFTXj+2XQmRwcrZ1bTAjLe8ZQHYoD8dQwGs+T+u+nYwo26U7mjobYL6wDRv2Xh/WMNV+Fl9/PxZdKqQjfN3nGN53OOaEFJWWQpLHc0/R1JgwHDsWjtwumbxdn2nY+/tcHElNh8n3PrzyTDO42N65s3K/3uW+VVAFvXcV9F6Z/fFYELZhE46q82iq2gU96hWNNhsFT1UiIqJ/E6M32o2cgmlTX0bvyurr07Eiur44C0u2HcGpkNM4FrQKP77VHzWvlqxcUL3Pq/hm+S4EnwrB6WO7sea/b6JfjRymAfCohXtGTcLsuSuwZc9hnFB/I+s9cWArlv4wGUOaeWX40jbCp8MoTJn2ET58+ym0zu5JlUtdPPjaVEz78B0MbpL34p6jX2sMffdHLA08iJOnQxESchR7N/6Jb995Aq18pfBihLfe9jS8P64Pcu5e64bGj03Eh9Pex0u9K2cqcOS+jXxwq4W+L8/E3HVBCD59GiGhpxC8ay1+/WgUOlbM52OkfJ4Dr/bP4IOPp+ODIU3gIk+kjdVx39sfYfr06Zj+8ceY8kRT5Drpw12brzzRevj7+Hj6VAxr6am3YSjfA+OmqbSZ/jE+mjQYjXXLcJWf2o1U25yKlyWfGH3QYewPWHfwJEJOByPwp+GokylLuKFW35cxc+46BAWfxumQUJwK3oW1v36EUR0r3nAgQdcaffHKl4uxaf9x9XencFTSdsL9qO2mkr32g3hz6jS881hj26dzo/a77XBM+vBDvDOsNa5LKqdK6DJ6OuaofHk0JAShp0/h8PaV+HnacLQpm026ZucW5w2PWvdg1KTZmLtiC/YcPqHWeUql+Qkc2LoUP0wegmaZDir/ecOz9XC8r66PqcNawlM+biiPHuOm6evj448mYbA1A1x3nEafDhj7wzocPKmOITgQPw2vk7mpe0Gvedca6PvKl1i8aT+On1bXwdEgrPpxAu63ZgA8+Kba93ceQ+Ncb5e34v7rgXYdm6KEIR1X9m5DUHbdUFxc4HBpJ/77cn/0fPIjrDyVQ1+VfLm5+5izf0+M/c8ibDl0Aqf1eViJ71/ri+rZHbI6z22HT8KHKk2Gtc543xFOqNRlNKbPUdf10RB1Tk/j1OHtWPnzNAxvUzbLZ+0MsJz8GZ//FQmzwQG1nngNg6tl/8nsGVGh+wt4T32XvT2kqbqzZMPoj94vv6/O41t4pH7mrixG73YYOWUapr7cG9bLsytenLUE247IdXQMQat+xFv9a14LtrhUR59Xv8HyXcHW63f3Gvz3zX7I6dadkVOlLnj+swXYuF/SWd2/967HnA+HouUNR+7M770yf/di87Ht2B1pBhxqoEPHKjmco9vL4OHhkW77mYiIiHLi2AwT1izG6GrRmPfCKBwbNBvj25dG4pkQRJp8UdXPQxW+zTj/11j0f24Dar72I754RhUI48IQcskZ5auUhQzKbj6/DOMGPIN5YRmftHnivq+2Yla/UqpwkA5zfCRCwi8ixdkHlSt7w9VkQHriAXz5+EP4IPCy/guX5q/jj3nPo76TBaG/DEWv19YhVr8jnNB43GIseKkxnKPX4PU+T+GXTNvLnlvjUfj6xzfR2dcBSLuMiNNnEGvyQiW9D0DKwZmo2vMTeD/yHTbOuAeeljD899FOeGtzNgXsUvfjmy2fo6+nGXs+6oF+nx+D7EHetvGxdR32NK8BBM/qj3um7snUPNtU6V5M+9+nGFRblQzNiYgKDcOl9FKo4F8W7pLWUZsw9YmnMXtvou0vbiS/58ARTd5YiT+fr5250pVB2rYpaP/Q17hh0t+t+cpYC2P/Wo1XG6vznI301M2Y2GoQfowwodmENVg8uhqi54/FyxdHYfbI+nCTCpZiuTgPw1qOxSrJYqZKuHfa//DpoNqqomVGYlQowi6lo1QFf5S1nnBsmvoEnp69F1nPuGfLcfjpp7FoqWrY6alxOBt6AfCpjPKeJiQc+AlvLq2M6a91hWHLRFR5+AfbX92Io22/deZE/3umYo89c6qK4eDv52FqD18YLSpfhoTgQqobfP0rokwJCw5+2ht9ph9SR5CLW5k3PO/DV1tnoV8pVR1JNyM+MgThF1PgrNKksrcrTKqCnXjgSzz+0AewZY185o3h2PHgH1j9amNkmwPSU7F5YisM+jEClgzHOX/sy7g4ajZG1nfT3ZfkyfDFecPQcuwqPTZLga95z5YY99NPGCtBNrXtuLOhuAB1HZT3hCnhAH56cykqT38NXQ1bMLH1QPxw/kYX7S24/zp1wIebfsHjFdKxf0Yv9P3kcDb5Q56cW2Cxr8pUD68sW46XGpgQ9+czaDb6LyTZ3sq7/N3HHJtNwJrFo1EDxzDn4w2o++xTaOyWhAuhZxDvUQlVvFzUXpoRsfI1PDBiDk5nPOwc7+dGVB/8PeZN7QFfdXyJUSEIuZAKN19/VCxTApaDn6J3n+k4ZE8QSavNKq3KW7B5QisM3tAfc5a9gzbuFkQsfg49X1iCqKvbVZX1J37FzqkdYIqah0qNx9qWCyd0+Wgr/jfYD+lBU9H1vlk4kfU0OffGf4K+wQOlzNgysTUe/uG87Y1raVEteh5eGHUMg2aPR/vSiTgTEgmTb1X4eeiLD3+N7Y/nNtTEaz9+gWeauiIuLASXnMujSlk36/W7bBwGPDMv0/dExnT+v/dWotYLo9GyVCqiw8MR7VAW/uU84SBBsKP/h2ceeROrL2TZ8QLdK+33tDzci4XRG0N+2oYPuzgjcfXLaDV0LqJtb90pRSG4QkREVHwY3NHhza8xvvE5/PxsD7Ro1QmdAtrg4Zk7cTndBL8+4/HxV//D7Geq4ti3w9G5ZRt07tgCbQd/jX0Jqrzh1wvjnm2bZQR7R7irGsilPXMxdcQ9aNagGTp064Fu7Zuhee83sDQ8DXCtj6FjH4R9dsSkXf/BOz8dR5oqGlUaOAEvtrr22MdUezjeGdkILulx2PLZFMzJrVAtPDtjwhcScDAh8eDPGN1VbbtTN3Rt3wSN2z2Kd37bj8vupdQHLYhatRj/RKt1msqjZ/822Tb/Ld2tDzp4SCVyL5b+dUIHNfK+jTxwaoQXZs9QFRwXJAX/jnG9W6BZ+67o2qEZmnYehR/2xcPg3R6vffoSWuRpzLj8noNU7Pt2NB4ZOAiPf7RBnXu1KG0vvhkxCIMGqdfAgXjo1d9wJg9Jr91t+cpyHD+PfQQDBw3Buyuj9Pk3n5iDlwdL+gzEwPtfwdyIzH9fosVYfDS8HkwhqzBz3AgMe34CPvtuNYJ17ccJjV6YjRmqoO6SFIzfx/VGi2bt0bVrBzRr2hmjftiHeIM32r/2KV7KesLdOuKNz8aghacBycfnYmzPFgjo2BEBzTrh6f9sx5W6T2LGyx1VahUOt54v4tXuvjAm7cGXg9qgWYfu6NG1DZo06YzH3/kBKw9EIV9PFW9F3nB0h5vpEvbMnYoR9zRDg2Yd0K1HN7Rv1hy931gKa9YYirEPVrxaWchf3kjF8Z/H6utjyLsrrRVN8wnMeXmwvj4GDrwfr8yNsN4XriqBFmM/wvB6JoSsmolxI4bh+Qmf4bvVwdYKcIGveTd0fOMzjGnhCUPyccwd2xMtAjqiY0AzdHr6P9h+pS6enPEyOuYjAxT2/dfoVx+1pbVaehJCjodkE9QQGYIahaaA9zFjDQx87SnUODMX4+9thebtu6BDs7Z4/Os9iFdp4ttjPMb28LB9OBduPfHiq93ha0zCni8HoY3Ki917dEWbJk3Q+fF38MPKA4jK4YIxGE0qW/0fPvq/Y0hT10LZe8fhuZbXzsXtYHDvgDe/Ho/G537Gsz1aoFWnTgho8zBm7ryMdJMf+oz/GF/9bzaeqXoM3w7vjJZtOqNji7YY/PU+JKi08us1Ds+2zeFLylgdj058Ds2S/8FHj3ZA8zad0LFlC/R6/S+EqWu0RK3H8MHr3eBp+7jVTdwrbW58L7axxOLY8fMqVxpQoqbKv4V1A70JDGwQERHlh6EEfEtH4PsRj+HtP49ZB0CzRGPb59MwP1wVR02V0KlXNRydORSPv7sCIfoRmgVRG2dgxl8RqsBqQnlVGG+YqRBwEXOfC0Cbfi9h1rKDiLpaeLAg7uAvePf7narQZoBr49ZoWsL2ltry1k8n47dQtU3HWnjy7WfQQDqyG6uoytPzaOkOxO/6QhW+bUGFGzKi8iPP4aEqDsCVPZj14kQsOXHtOU5i2EZ8N64vWnd/z7ogZi0Wrrmgj8W3ez+0u64Nrxe69W4Pd2M6kncuwRL92C6f27ghI7z7vYDhjdWGE7bjk1Gv4veDsVePM/HkUkwa8yX2JqukqfkInuiabSPjLPJ/DixRwQjctBGbj0bBLAVvVZE5tX0zNm7ciI2bNmHHiWv7lKu7Ll9ZcPHoNmzauBlHLqTqinx6Qjh2S9ps3ITNu0Otx3iVEe5VqsInchlef2w4PvptGVYs+hGffL4UMoyA0bsfXhjeWFVRE7D9k1F49feDiL12wrF00hh8aT3heOSJrhmalRvhO2AEHvBXlca04/jvq29h/lHblpNCsOrDJ/Dk9O1INMnsEoXBEbWbN0YZVcJO27MA3wdagzrCcvkE1n83GTNWROaSdlncirxxcS6eC2iDfi/NwrKDUdbAgbDE4eAv7+L7nWlIN7iiceumuJo18pk3LBePYptcH0cuINWaARC+W87/RmzavBuhmTOAWo87qlT1QeSy1/HY8I/w27IVWPTjJ/h8aYhaX8GveaPvAIx4wF+lQhqO//dVvDX/qC3vJSFk1Yd44snp2J5ogkO+MkBh3n8Bh+o14C/bt0Qg7DZPFVOg+5ja1+TD32PEwFfx675o6/vmSPw9/QPMD5M86Y2OPVrqj+bGsXZzNLZeMFjwfeC11haWyzix/jtMnrECkTkloj5nSdj55XT8EWGGwaEGBr/2BGrexhquoYQvSkeotHjsbfx5zJqzLNHb8Pm0+bBenp3Qq9pRzBz6ON5dEWJtVWOJwsYZM/CXDKZrKo/2XRvqv7uOQR1g1BpMeHQYPt90Fil6YQIO//o63v0zUl/b5e59HPdmmDu34PdKuxvfi68xIyzsrD73xnLVUDUvX7O3GAMbRERE+WJG6Jy3MW1TTObCXvIe7NiXqCtwaYe+xxuf7YStBbdNAoKCjuiCo7GcP/ylEJxBalwcErMtvFkQefwUotXfGVy84KM7q9vE/oOP31uEc2YDXBqNxDtP1kKlhyfipQ4lYUjaj2/f+RZH8lJGNpZHz3tawEWatu6Yh7lHs/sjCxIT7Q2dE7Bx4UqclTJZ2e7o3zFzicbo3RV92nnAmJ6IbUuWWZ/25XsbN1ISXfp2RCljOuI3/orfjqkdycJ8fBlWHVbLjaXRNKC+bemNFegcFJq7MF/ll/kcFr39JhaEXH8+S3bpi46ljEiP34hffzumUisL83EsWyXN940o3TQA9a9W4j3QvkuAbk5tPrQYv+3K2kklAXtmjcH0rdY0vnkWpKRYAznGCrVQt1AK+7cib6QiLi4x8/rsLJE4fiparVedfy8f6/gYdrc4b5jPLcLbby7A9Vmg4Ne8R/suCLBmACz+bdd13ZQS9szCmOlbkZjfDFCIaeHi56uOTf2QfglRF64/tiLHHIpfJ07FhotZclDSbmzdnaDyjhFl/KvYFt6YJSXFGvgyVkCtAl4wlovLMeOLrbhsMcCt5Wi8cp/P7avkqrSY8/Y0bIrJnBbJe3Zgn85UaTj0/Rv4bGfmqxMJQQg6YtZpVc7f37YwC0sUlk1+Bb8cs4Y0rrJEY92idZCsYnBthjYtr7W8KPi9MoMb3IuvsSA6KhrJ6hANJm/45Xd8rFvgtp1zIiKiu8WVpCvZVAjSEBsTb12enJRtITkx9rJ+4mJwcIWry42+go1wLlUelavXQp26dVHHTxX2ZH0GExwzFUIsiFrxAaYuj4TZ6IHWL/+AX9/sAW9DCoJ/fNv6VCYvHOuhQS15Ym1B2L59yFpWzU5S4EIsk0KPsQw69+uSoSmsKtB264u2HgZVsNqKJculqapSgG3kyKkumtRz0eu6EHUFlevVQ72srzpqj/RJMMKvYgXr3+VLXs9B4bnr8lU+mU/9gZ9WX2vhcI0T6jappwc2tFyIwpXK2ZzvenVUHrRW8o1+FVHBfjyONVCvpuQVM6IO7MXp7MrpqiJ/NiJF/+3NM+Po3xsRkpYOU+XH8PkvH2BIgJ86gptzy/OG0RmlyldG9Vp1ULduHVzLGo5Zuujcyrxhxqk/fsLqawMkXFPga94RNerV1HnHHHUAe7PPAIg8G4GUfGeAwkoLI0qW9IBJ7SMsiYhPyOb4C8wRpSurc3pdetVFrfIFCyJYXUHSlez2Mw2Xoq0tPAyu7tZFuTAf/RsbQ9KQbqqMxz7/BR8MCYBfvi8YC07/Mg3/DU5FutELvV5+Aa1v5vDyJYe0SItFTLx1eXJSdoHTRMRe1lcnHFxz6D6Tfgkhp7Nv+Zd84ACOyxsq71WpXt668GbulRnkfC/OzJJwGUn6RuEOz1sS8M8fBjaIiIgKiTk9l5Kxet9aUDCo/7JShdt6A/DK5/Owfncwjh3Yga0b1mPtmjVY9nF/lM3pYYglAovf+xjrL1lg9KyKat5GpJ38P7zz+Y7rnkzmyLWs7Ym9BRcjpXlrHqQGYeGSYKSlG1G6U390sw+NoQqVPfq0gbvBgrhNf2KVfVCzgmwjJyY/+OrR4E2oNuRrLF+9Gquvey3Ba22tJTYH5xz6L1+ngOfgFiu2+Sq/zOYc8oUJfr7eutBqqjYEXy/P7nyvxpLX2lor4Q7OcLaPWGn0gY+XpEo6oiKkW8atl7JjJsb/ZyeiLUZ4NX8S0xZuwbblX+O1h5rCu5Dz0E3lDWNJ1BvwCj6ftx67g4/hwI6t2LB+LdasWYaP+5dVqZ6DW5g3zCoPZKvA17wRPj5e+tjToyIgLf8LVSGlhbOTvXaZhtQsD+dvilt/fKyu9TXXpZVa9ts4tMimUntzVH4z2+75eZ2WNGUHZo7/D3ZGqzT0ao4npy3Elm3L8fVrD6Fpfi6Y5CB89bG1BY1D1Ufx2lO1bHk4l2vkljHL5XdD6em2tJIuJ/lkuXwBF61RBXiWsn8B38S9MqMc78WZpaelqTKA/OQIJ8f8H0Nhy2OOIyIiolvHiEoDZuCvJbPw0oNtUbN0CsL2bcH61SuwfPlyrAgMsT4VyYEl7jwi4q99wJKShCtptl/ywOjkBCdbmSQtz6VqMw4tWoL9qlRjLNkB/XpYC1NG757o09oNBksM/v5jDS5aP1zAbeTA6ABH6+NNxARvwAqVRpJO2b6WzsMXP2+y/t0N3dw5KJrubL4qPEY4OJp05dQSE4wNK7I5z1dfSzHvi5+xKd76l5JXHHReUX+bU6W50MUicPojuOfRSfh5Uyji051QttG9ePGzP/DPknfRxz5S651krIQBM/7Cklkv4cG2NVE6JQz7tqzHap22KxAYknTD6uBtzxsFvuZV3nGw5h2VAW5JYKsw0iIl1d5nxRGON9u8J6Mr6nr5c1m26bRsxU7dnbAoiA2cjkfueRST1HkLVWnpVLYR7n3xM/zxzxK82+faALa5iVk9A59viIPF4IpmI8ahv4/tjbuRJQ2p1qgCTEZ7AOgm7pUFYHBQ91d9caUixbYvdxIDG0RERHeY0fs+TJjyIKqVSEd04Ew83q4p2vV5GEOGDsPw4cMx+rMNkElIsueB9q9MwsBKRpijzuOCKlA71Xkak0c3ynPzd0tSAhItUigxwM0j8/jqN2I5sRiLdyYj3eCBdv16wttohE/P3mjlZoDl4jr8ue7aBIgF3Ua2LAlIuCLrsuD8yg8xUqWRpFO2r5EvYfrqcOvf3cDNnYOi6U7nq8JjQULCFV3RtpxfiQ9HZnOer75G4qXpqxFuPy5zIhKT5BcD3D09b2PBNwVnNn+H1we2Q+ueI/Hhgv2IsRhRqvFTmDHjKVS5oyVwI7zvm4ApD1ZDifRoBM58HO2atkOfh4dg6DBJw9H4bINtQMhs3YG8UeBr3qzH7dE5wN0z83ghhaIw0sKC+MsJ1vQ2usLDvRB30nwAP70yMtt0euaDZTib80m+7VLObMZ3rw9Eu9Y9MfLDBdgfY4GxVGM8NWMGnsrrBWMJw+8ff4+DyekwevXCC6Oaw0l970jOyc61FhXW4Gex4uQBD+lzoo7ucpz9u/Ym7pUFYHTz0N1eZGDg+Lg7n5nu6G2ViIiIAOeAbuhQWn0lpx3FT5M/wfozeW/R4NbqJUx5vBocLBH4462H8crccFWUd0HjUZMwvHYem/FeCccZ3WXEiErVqmWZMvQGLOFYsngrEtINcGt9L3qV80PP3gFwNZgRueYPbMg460FBt5GdlDCc1o8aZdC1ShlmbSi4mzkHUrAsiu54vio0KQg7fVY/bZdBMCvl54SnhSP8vJwfI8pW9kf2PdlNMJluVZHYgujDy/CfFx/AYzP3IkldK+4BfdC9/J0sgjsjoFsHWLPGT5j8yXrkJ7vfkbxR4Gs+DeHh5/UVaixbGf45DGVgMpnUmvOvsNIi/kIULsvt0VAGPj63+/qyKxr3MUv0YSz7z4t44LGZ2JuUDoN7APp0t48hkbuUvd/gowXqPBgcUHPwWNznk1OFOx1Xkm3joJQoYa2gZ+XgBCfjnQx6GJHT5k3Va6KKNINMT0LYKXsg7ybulflmRGmf0nCWXbBcREQRGPT2Tt5ViYiISHHwsD/1uIDzZ68vXLqVKaULD9dxaYbnpwxFLcd0XFr3CaYtO4l1H3+IZaqAYXBviRcmP4FqefmmTzmAHXvliaERZdp2R8tsC/9GuJUqafvZzoILyxdhY5wFBtdW6D94CO4NcIXBfBarFm/O3Me8wNvIhjkYgTtlnA4j3Fv3QIc8/EluCnwORGqqdbpMgxNcXPSSIuGO56tCY0Zw4E5EqnKz0b01euTnhFtCsPeATE1sgGvzjmjrYVuegVvL5/Fcdw+Vm26lRBz6J1A/ITcY3FDSI6fMdDs4wMNDBuKUrHEe12cNN5Qp5Zz9M+w7lTcKfM1bELL3gG32iObomH0GwPPPdYdHfve7ENMi9dQp65Nzoy8qVrr9baK0InYfSzz0DwKtFwzcSmZz3nIUh78/m4l1MmaHZweMGFwvh/FizLgQedEa9KpQHdWu+07yRNtXXkSPO3mtGqvhgdH9Uem6fGRCne5dUF0tT0/ej2077bOJ3cS9Mt9MqFypvL5vWiJDcDrLpC93Qn4vYSIiIipkyWFhiJBCraku2nf0yvDlbETpgDH4ekofyDT/mTmh0eh3MaK+MxC/HbM++F0XjC2Rf+KDTzci1mKEZ7uX8PbASnn4so/F2gUrdGHIVHEAXn22Wea57Y2eaPTEF1ix9l3bgmssF1dj0bpLsBicEfDsaLRWNWlz6Aos3p51RoCCb+N6ydj66zwEp6bDVPY+vDmxDypkW3J1QYUGdeCXhwegBTsHVmmREdATORiro37D2zYUf67ufL4qPMlbf8U8mfHAVBb3vTkRfbI/4XCp0AB1Mp3wZGxbskaPJWD07oMXRjfPlO/cGg3Dl7OfQ0OZDrQwGCug/xvv4ukA65gz1ziiaqvmkIYaloQQHA9XO3THJCMsLEJV+SVrtEdHrwx7aiyNgDFfY0qfMtmc3zuZNwp+zSdvW4I11gyAPi+MRvPMGQDDvpyN59R1m78cULhpYQk7guPSlN9QAtVqVs2hIn5r3Zn7mBEV+r+Bd58OgB4bNgPHqq3Q3HrBIOR47t0JM7KcmY+PvtuHK+km+JbzzSE9LQjbtdd6b/DohCEDa1gH1NTc0GD4F/hiWC045jYC6K1kMMGv9wf4bnIP+GZIH8fqg/H6E3XhYLAg9p/5+CtDn6KC3yvzyVgSNar7qTOYjrSTRxCcz3FlboXCv+8QERFRvqTsnI+Fh5ORrgrevVUB+Zu3n8NTTz2D8TPmYuWccWgUuQuHsvRfdaw7ElNGNVbF+GQc+G4K/nvUXlGyIHTOu/gyKEFPe9f11bdwX8YSUQ7iVs/AtFUya4Q7mo2Zg2X/ew8vDR+Kp599AzMXrMPiD/rBP+XamBnXxGH9wtWIMBtgcnRSBS0zTi7/A0HZNG0v+Daul7p3NiZ+vR/xqihafdBXWLH8O7z78kg8OWQIho4YgzenfYuFm4KwZelnGFwr98JbQc6BnTl4EwJ16dgL/SZ+gQkjh2L4y+9j1hu9dHP/O6Uo5KtCk7oXsyd+jf3xUqgfhK9WLMd3776MkU8OwZChIzDmzWn4duEmBG1Zis8G22dDsErc/CU+XxcFi8EFjZ7/CQu/fAujhw3HmPd/wop5k9DV5QC2HVbHZfv8TTF4oU63JzBl/kas+/kjjB89FE8MHY3XZs7HnFdb6HQN/r/vsfqOPt1Mwc75C3FYxiHw7o0Pfv8Gbz/3FJ56ZjxmzF2JOeMaIXLXIWTN7nc6bxT4mk/cjC8/X6cq7Qa4NHoePy38Em+NHobhY97HTyvmYVJXFxzYdhgJ+cgAhZ4WybsRdEhVRlXOrdKgIexzXGRk9OqD95esxbp166yv1d/g8RrW43Tr/DaW2Zev/ROTe+b/Sf2duY8Z4FWnG56YMh8b1/2Mj8aPxtAnhmL0azMxf86raOGikib4//B9vi+YVBz69gP8ciotw3V9/QlO2f5/+GW3nCtPtH1rDn6d9hKGD38Rk39chgXvdAL++RTfb79zNXbLpV1YtycN9Z76Dqv+/BxvPDsco179BPPmvYsuXkb1/mZ8Om2hNYBtdxP3ynxxqIeGtR3VGVTf+buDcCn7r8fb6g5+3RIRERUjlgTExiTBbE5CTGyi+irPKh2J0TFIMZuREB2LK7alGZnjoxGbakZaXDTiUjOUAlL24LMXXse8w3FId6+N3qPexHvvTcSLD9VH/LK3MejRj7DlohnmpHjE64YQHug+6mk0dbUgbv+3mPTlXlW8ziD1CL6fPBu7YtJg9OmOYQNr2d64AUsY5o4Zgtfn7kd0uhtqdH8Kr0x+H+++9TweaumF6F2/4PURH9k+nFniprlYfFIK5SoV0oKx9I991ibNWeV3G5ZExNnTPC5RVRkyuozAaUPw+ORFOBidjtL1e+Ppce/gg2nT8P6k8XhuSB8EVErDibVLsDH0+rN1nXyfgwySAzFryhwEJwBOlbth9DvvY/K4J9G7sW/uAwje7flK7X/cpWi9/4mXYrKZhcWChFi1XL2fFBOLxBucqsuB0zDk8clYdDAa6aXro/fT4/DOB9Mw7f1JGP/cEPQJqIS0E2uxZGNo5nS0hOK3sSPw4ZowJBtLocF9z2LClMkY/2QXeIcuwBuPPonvD1vzb97lsN/mA5gz40cERjqhRtfBGDPhfUx9f4JK8xbwtURg67fP4+lp2/M2FegtzBspez7DC6/Pw+G4dLjX7o1Rb76H9ya+iIfqx2PZ24Pw6EdbcFGOLT7elgcKnjfS4y4hOkXls8RLiMluGp5cj9OuoNe8BaG/jcWID9cgLNmIUg3uw7MTpmDy+CfRxTsUC954FE9+fxipec4At+L+ewFbNh+GOd0A58at0Dy7riAlfFCtdm3UvvqqCu8S0s7EAJNnBdTK8F5V7wJ0Z8nDfcySGKfOodwzYhCb7cVqQWJMDJLleld5MpNsz7MZB+bMwI+BkXCq0RWDx0zA+1Pfx4QXH0ILXwsitn6L55+ehu0ZL5j0BFyKToE57TKi1b85StiM6a9+hs1n1LbS0xB96IDtjQzMRzD7hXH4v/2xsDiWR9shr2Dy5NcwrJsXTvz6EgY+8wP2X0xQ21L7Hpf5W82SEHvjtEhPRHSM2k9zAqJjs706Ea/SKNWchrjoONsyK0nn2CtmpF4IxOdDH8NbC0+ofPEAnn9rMt4eOwgtyxqRcGIppgwdhR+Cr/+2Ldi9Mu/3YmGqEYCmMl+4+RwCNx29wXV7+xg8PDzydx8nIiKiW0NVumq27YzW9crDLTUSR7euxd9HbjQ7wa3hWqEpOrRrDH9vFz2V4Ymgjdh0KCr7YIXmhQe/3YDP+pSEZd8M3HPvDBzJpZST/23cgIsfGrdpi4bVy6GMiwHJl6Nw9vQR7A3aj9DYfBa3buIcOFdojm6dm6JKGSekXDiOwLXrceBCgY6ocBWRfFV4XODXuA3aNqyOcmVcYEi+jKizp3FkbxD2h8beoIDtjPJNu6Jz82oobYzHuSNbsHbTMcRanNH7P0H45oFSMG+ZiCoP/2D7/E1w9EKtVm3RvGZ5lSctiD13HEGbNuFQUcgPGRhL1UTbzq1Rr7wbUiOPYuvav3GkOEz/U8Br3rl8U3Tt3BzVShsRf+4ItqzdhGOxFjj3/g+CvnkApcxbMLH1QPxw/vangan+K1i+bCzqG8Lx46COmLAlPwMYF547cx9zhFetVmjbvCbKq2vaEnsOx4M2qe+ECwX7Tsgvkzfqde6C1rV84Zh4Fgc3rcaWE5eL0D3ShNJ1u6BXu9oobYjD+WPbsH7jUUTn+vVW0Htlbozwf2YB/p4QANP5X/Bkh/FYbx/m4w5iYIOIiIhuitF3CP634UN0dTNjz8c90G/m0WJcaaZ/n1sQ2KBipSgENmCqhTF/rsKrTYw49c1D6Dp5++2p1BPll7EShv/+Nya1dUT4/x5Dpzc3ZW61dIewKwoRERHdBCMq9huAtm4GpKftw4qlxxnUICLKL/NxLJi7BQnpJvj3vg8BNzUnNtGtY6zSD/e2UBk07SgW/ra1SAQ1BAMbREREVHDGGhhwfws4G9KRunc5lp5kWIOIKP8sCF/0I/46Z4apQl8M7OZpW05UlDii4SMPoqkjEL/5J8w5UPAOLYWNgQ0iIiIqMFO9AejXwAGG9BQELVuK04xrUHHEjtn/ckUkA8StxeyvNuBCSkm069batpCoCDHVwz1d/GG5fBA//Wc+worQdz7H2CAiIqICM5YNwEMDWqBU4hGsW7QOxxNsbxAVG0Z4NxuAAa18cOXAMvzfxjDbcvq3MHo3w4ABreBz5QCWzdmIsKLzEJqI8oiBDSIiIiIiIiIqttgVhYiIiIiIiIiKLQY2iIiIiIiIiKjYYmCDiIiIiIiIiIotBjaIiIiIiIiIqNhiYIOIiIiIiIiIii0GNoiIiIiIiIio2GJgg4iIiIiIiIiKLQY2iIiIiIiIiKjYYmCDiIiIiIiIiIotBjaIiIiIiIiIqNhiYIOIiIiIiIiIiing/wF9xs5c2ODoPAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnorm = matplotlib.colors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)  # min, center & max ERDS\n",
    "for subject in subject_list:\n",
    "    for condition in epoch_event_dict.keys():\n",
    "        con_epochs_path = f'{epochs_path}/{condition}_{project_name}'\n",
    "        # Get file for this subject\n",
    "        file = [x for x in {con_epochs_path} if x.startswith(subject) and x.endswith(f'{condition}-epo.fif')][0]\n",
    "        # Load epochs\n",
    "        epochs = mne.read_epochs(os.path.join(con_epochs_path, file), preload=True)\n",
    "        \n",
    "        # Calculate ERSP and ITC\n",
    "        power, itc = mne.time_frequency.tfr_morlet(\n",
    "            epochs, \n",
    "            freqs=np.logspace(*np.log10(freqs), num=(freqs[1]-freqs[0])), \n",
    "            n_cycles=tfr_morlet_eeglab_n_cycles(morlet_win, freqs), \n",
    "            return_itc=True, \n",
    "            average=False\n",
    "        )\n",
    "        \n",
    "        power.crop(wtimes[0], wtimes[1]).apply_baseline(baseline)\n",
    "        power.save(f'{ersp_path}/{subject}_{condition}_power-tfr.h5', overwrite=True)\n",
    "\n",
    "        itc.apply_baseline(baseline)\n",
    "        itc.save(f'{itc_path}/{subject}_{condition}_itc-tfr.h5', overwrite=True)\n",
    "        \n",
    "\n",
    "\n",
    "        \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mne",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}