Fig2_C_right_G_middle.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d831d9eb",
   "metadata": {},
   "source": [
    "imports and upload data\n",
    "-----------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ce2b0e40-1047-444a-90cf-11299377cceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "from tqdm.auto import tqdm\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n",
    "warnings.simplefilter('ignore')\n",
    "import os\n",
    "import pickle\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score\n",
    "from matplotlib.gridspec import GridSpec\n",
    "from collections import defaultdict\n",
    "from sklearn.metrics import confusion_matrix\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "07fbcce1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# upload data\n",
    "cpath = os.getcwd() #assumes mean_df_females is in the same path as the notebook\n",
    "mean_df = pd.read_csv(cpath+'/mean_df_female.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f17ff0d",
   "metadata": {},
   "source": [
    "# ----------------------------------------------------------------------- \n",
    "# ANALYSIS\n",
    "# -----------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1be321b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot/colors definitions\n",
    "data_loc=cpath+'/figs/'\n",
    "try:\n",
    "    os.mkdir(data_loc)\n",
    "except FileExistsError:\n",
    "    # directory already exists\n",
    "    pass\n",
    "\n",
    "plt.style.use(\"seaborn-paper\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c28bf6ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _plot_cm(y_true, y_pred, ax, ax_labels, title):\n",
    "    cm = confusion_matrix(y_true, y_pred)\n",
    "    cm = cm / cm.sum(axis=1, keepdims=True)\n",
    "    im = ax.imshow(cm, cmap='copper', vmin=0, vmax=1)\n",
    "    plt.xticks(range(len(ax_labels)), ax_labels)\n",
    "    plt.yticks(range(len(ax_labels)), ax_labels)\n",
    "    ax.set_xlabel('Predicted')\n",
    "    ax.set_ylabel('Real')\n",
    "    ax.set_title(title)\n",
    "    return im\n",
    "\n",
    "def plot_cm(y_true, y_pred, y_shuffle_true, y_shuffle_pred,name):\n",
    "    '''\n",
    "    plot confusion matrix\n",
    "\n",
    "    Args:\n",
    "        y_true ([np.array]): array for true label\n",
    "        y_pred ([np.array]): array for predicted label\n",
    "        y_shuffle_true ([np.array]): array for shffuled label\n",
    "        y_shuffle_pred ([np.array]): array for shuffled predicted label\n",
    "    '''\n",
    "    fig = plt.figure(figsize=(23, 10), facecolor='white')\n",
    "    gs = GridSpec(ncols=3, nrows=1, wspace=0.1, figure = fig, width_ratios=[10,10,0.3])\n",
    "    fig_ax = fig.add_subplot(gs[0,0])\n",
    "    labels = np.unique(y_true)\n",
    "    _plot_cm(y_true, y_pred, fig_ax, labels, f'Real Accuracy {accuracy_score(y_true, y_pred):0.2f}')\n",
    "\n",
    "    fig_ax = fig.add_subplot(gs[0,1])\n",
    "    im = _plot_cm(y_shuffle_true, y_shuffle_pred, fig_ax, labels, f'Shuffle Accuracy {accuracy_score(y_shuffle_true, y_shuffle_pred):0.2f}')\n",
    "    fig_ax.set_ylabel('')\n",
    "    fig_ax.set_yticklabels([])\n",
    "\n",
    "    # plot colorbar\n",
    "    cb = fig.add_subplot(gs[0,2])\n",
    "    fig.colorbar(mappable=im, cax=cb, label='Fraction of labels', )\n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "    fig.savefig(data_loc +name+'.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "17ec3619",
   "metadata": {},
   "outputs": [],
   "source": [
    "def just_average(df):\n",
    "    d = df.mean()\n",
    "    d = d / d.sum()\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "530ebb25",
   "metadata": {},
   "source": [
    "# Phase decoder avaraged per phase per mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "efe893a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# decoder parameters\n",
    "df=mean_df.copy()\n",
    "it=1000\n",
    "leave = 'mouse'\n",
    "predict='group'\n",
    "minum=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f7b794ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose mice that have at least 3 phases and 10 repetitions each\n",
    "minp = 3\n",
    "minr=10\n",
    "phase_counts = df.groupby('mouse')['group'].nunique()\n",
    "rep_counts = df.groupby('mouse')['rep'].nunique()\n",
    "keep_mice = list(set(phase_counts[phase_counts >= minp].index) & set(rep_counts[rep_counts >= minr].index))\n",
    "\n",
    "# use avarage per phase per mouse\n",
    "data = df[df['mouse'].isin(keep_mice)].pivot_table(index=['group', 'mouse', 'uuid'], columns='syllable', values='usage').fillna(0)\n",
    "data_avg = data.groupby(['mouse', 'group']).apply(just_average).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2c3931a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b1533d959c62447eb560c3e0630ecf19",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e5f30afc6644e29821d07a51ebfad1f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#decoder\n",
    "clf = RandomForestClassifier(n_estimators = 250,class_weight='balanced')\n",
    "out = defaultdict(list)\n",
    "for i in tqdm(range(it)):\n",
    "    sm = data_avg.copy()\n",
    "    a,b= [],[]\n",
    "    for leaveout in sm[leave].unique():\n",
    "        train = sm[sm[leave] != leaveout]\n",
    "        test = sm[sm[leave] == leaveout]\n",
    "\n",
    "        X_test = test.drop(['group','mouse'], axis=1).to_numpy()\n",
    "        y_test = test[predict].to_numpy()\n",
    "\n",
    "        #balance labels\n",
    "        btrain = train.groupby(predict).apply(lambda x: x.sample(n=minum-1,replace=False)).reset_index(drop=True) # new balanced df\n",
    "        X_train = btrain.drop(['group','mouse'], axis=1).to_numpy()\n",
    "        y_train = btrain[predict].to_numpy()\n",
    "\n",
    "        clf.fit(X_train, y_train)\n",
    "        y_hat = clf.predict(X_test)\n",
    "        a.append(y_test)\n",
    "        b.append(y_hat)\n",
    "    out['y_true'].extend(a)\n",
    "    out['y_pred'].extend(b)\n",
    "    out['acc'].append(accuracy_score(np.concatenate(a),np.concatenate(b)))\n",
    "\n",
    "#shuffle\n",
    "for i in tqdm(range(it)):\n",
    "    sm = data_avg.copy()\n",
    "    sm['group'] = sm.groupby('mouse')['group'].transform(np.random.permutation)\n",
    "    a,b= [],[]\n",
    "    for leaveout in sm[leave].unique():\n",
    "        train = sm[sm[leave] != leaveout]\n",
    "        test = sm[sm[leave] == leaveout]\n",
    "\n",
    "        X_test = test.drop(['group','mouse'], axis=1).to_numpy()\n",
    "        y_test = test[predict].to_numpy()\n",
    "\n",
    "        #balance labels\n",
    "        btrain = train.groupby(predict).apply(lambda x: x.sample(n=minum-1,replace=False)).reset_index(drop=True) # new balanced df\n",
    "        X_train = btrain.drop(['group','mouse'], axis=1).to_numpy()\n",
    "        y_train = btrain[predict].to_numpy()\n",
    "\n",
    "        clf.fit(X_train, y_train)\n",
    "        y_hat = clf.predict(X_test)\n",
    "        a.append(y_test)\n",
    "        b.append(y_hat)\n",
    "    out['sy_true'].extend(a)\n",
    "    out['sy_pred'].extend(b)\n",
    "    out['sacc'].append(accuracy_score(np.concatenate(a),np.concatenate(b)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2db872de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save decoder results\n",
    "decoder ='phase_decoder_averages'\n",
    "with open(data_loc+decoder+'.dat', 'wb') as f:\n",
    "    pickle.dump((out), f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3672f791",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABU4AAAJbCAYAAAA7eQ5eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzde7RWdZ0/8Pe5gIAJ5iVvJRikyUUOdFBAUSHIy1B5q4xR1C402DgOpTUuqIyWqdXKdNZYeSmxqQzL0NKxIkVRIUVEBFKXF3TMO4qIgoBn//5wPD95RNxwOJwDz+vlelY8z97f7/48RxZtP7y/311TFEURAAAAAACa1bZ1AQAAAAAA7Y3GKQAAAABABY1TAAAAAIAKGqcAAAAAABU0TgEAAAAAKmicAgAAAABU0DgFNpsOHTqkoaEh++67b8aMGZPVq1dv1Dxnn312fvKTn6zz2NKlS9OxY8dceeWVLSl1syuKIv/yL/+SXr16pbGxMQ8//PDbzrnzzjvT0NCQhoaGNDY2ZtasWc3HzjzzzPTt2ze9e/fOhRdeuDlLBwCqxOTJk9OnT5/069cvjY2NefTRR5Mku+666wbNc8stt6RPnz4ZNmxYVq5cmREjRmTAgAG5+eabN3iuJPnUpz6VESNGbPC4tvbHP/4x++yzTz70oQ/lsssue9vxJUuW5KMf/Wh69+6dfv365eqrr17reFNTUw444IAcf/zxm6tkgKpT39YFANVjxx13zLx58/L6669n1KhR+c1vfpMTTjhhk15j2rRpaWxszNSpUzN27NhNOvebXn/99dTV1W3SOa+//vosXbo0Dz30UK6//vp8/etfz29/+9u1zunbt2/uvvvu1NXVZcGCBRk7dmzmzp2befPmZfbs2Zk/f35WrlyZffbZJ6ecckq6du26SWsEAKrXHXfckVtuuSXz5s1Lhw4d8sQTT2TbbbfdqLl+9atf5ZxzzslRRx2VWbNmpWvXrrnppps2aq5XX301d955Z7p06ZJnnnkmu+yyy0bNsz6tce+3Zs2anHHGGZkxY0a22267NDY25uijj86OO+7YfE5tbW3OO++8DBo0KM8++2wGDhyY0aNHp3PnzkmSyy+/PD169NikdQGwNolTYLOrq6vL/vvvn3/84x9J3rhxPP300zNo0KA0NDTkD3/4Q5Jk1qxZGTJkSAYMGJARI0bkqaeeete5p06dmu9///u5//7789JLLzV/fumll6Zfv37p379/vvWtbyV54z8A9t9///Tv3z9HHHFEkuTkk0/OjTfemCRZvHhxBg8enOSNlOspp5ySIUOG5Gtf+9o71rZs2bKMGTMm++23XxoaGjJnzpwcf/zxa/3HQGNjY/N3f9N1112XE088MUly5JFHZvbs2SmKYq1zunTp0nzT/sorr6Smpqb52MqVK7N69eqsWLEi2223XbbZZpt3/VkBAJT19NNP533ve186dOiQJHn/+9+f9773vc3HTz/99PTp0yef/OQn8/rrrydJevTokZUrVyZJrrjiivzHf/xHrrzyylx99dX56le/mi9/+cs58cQTc8stt6ShoSFLlixZ65rnnHNOBg0alP322y+XXnrpOuv64x//mMMOOyzHHntsrrnmmubP77///hxyyCHp379/DjjggLz22mvrvE+bMWPGWonNQw89NPfff38WL16choaGfPrTn07v3r2TJKNHj85HPvKR9OvXb61rVd5n3njjjWv9Bf6FF16Yc889d62677zzzvTr1y+77bZb3vOe9+Sf/umf8uc//3mtc9773vdm0KBBSZL3ve992XHHHfPiiy8mSV544YVcddVVGTdu3Dv+OwOg5SROgc1u5cqVmTVrVvOS8ssuuyx77bVXLrzwwixbtixDhgzJ4Ycfnj59+uS2225LXV1dpk6dmu9973u54IIL3nHeF198MQ888ECGDh2a0aNHZ9q0aTnppJMyf/78/Nd//Vduv/32dO3aNS+88EJee+21jB07Nn/4wx+y77775oUXXnjXuh9++OHceuut6dChQ5YtW7bO2iZPnpy99947v/rVr7JmzZq8+uqrOfnkk/OLX/wiI0aMyIIFC7LDDjtkjz32WGvuJ598svmzmpqa7LjjjlmyZEl22mmntc6bOXNmvvSlL+Xpp5/ODTfckCRpaGjIIYcckt133z2rVq3Kj3/8Y41TAGCTGjVqVL71rW+lT58+GTlyZE488cQ0NjYmSZ555pkce+yxufDCC/Pxj388N910U0aNGrXOecaOHZubbropxx9/fA4//PDMmDEjP/nJT3LVVVetdd6NN96Y5557LnfddVdWrVqVYcOGZfTo0dltt93WOu/qq6/OuHHjsvPOO2fChAkZP358kuTEE0/M9773vQwfPjwvvfRSOnTokIkTJ77tPm3u3Lnv+J0XLlyYX/7yl+nTp0+S5Morr8wOO+zQfL969NFH57777nvbfWa3bt1y+umn55VXXsm2226bX/7yl/nd73631txvvfdL3mhEV/7F+lvdc889ef3117P77rsnSSZOnJhvfOMb73g+AJuGximw2SxZsiQNDQ157LHHcuihh6ahoSFJ8pe//CWLFi3KFVdckSR5+eWX8+STTyZJTjjhhDz66KNZs2ZN9txzz/XOP23atHziE59ITU1NjjvuuJx33nk56aSTmpMEby5d32GHHTJ//vz06tUr++67b/Nn7+aoo45qTlm8+OKL66ztpptuam5o1tfXp2vXrhk1alQmTJiQFStWZMqUKc3J0o0xbNiwLFq0KHfddVcmTZqUP/3pT3nooYfy+OOP58knn8zSpUtzyCGH5IgjjlhrqRcAQEtst912ueeeezJjxoxMnz49o0aNytSpUzNq1Khsv/32Ofjgg5MkAwYMyOLFi1t8vb/85S/5wx/+kBkzZiRJXnrppTz88MNrNU5feeWV3Hnnnfn1r3+d+vr6PPnkk3nmmWfSuXPnvPzyyxk+fHiSpFu3bknWfZ+2Pvvuu29z0zRJLrjgglx33XVJ3liZ9PTTT6/zPjNJjjnmmPz+97/PgAEDst122+UDH/jARv8sli5dmrFjx+aSSy5J8kYT9cUXX8yhhx7a/PMBoHVYqg9sNm/ucfroo4/mwQcfzLRp05K88WCkyy+/PPPmzcu8efPy+OOPp3v37vnmN7+ZY445Jvfdd19+8Ytf5LXXXlvv/FdffXWuuuqq9OjRIyeccEJuvvnmLF26dINqrKurS1NTU5K87XpdunRp/vWG1FZXV5ejjjoqv/vd73LttdfmmGOOeds5u+++e3PKoCiKLFmyZL2Nz0GDBuWpp57Kc889l9///vcZMmRIttlmm+yyyy7p379/Fi5cuEHfGwDg3dTX12fkyJE577zzMnHixFx77bVJstZKl9ra2ual+uu7r3o3RVHkO9/5TvP94aOPPpqDDjporXOuv/76vPDCC+nVq1d69OiRZ5555m3Jznfz1hor63zrvd/NN9+cO++8M3feeWfuvffe7LXXXuv9Tm+uOLryyivX+Zfmb733S5InnniiOU36VqtXr86xxx6bCRMmZOjQoUmS2bNnZ+bMmenRo0eOP/74/OEPf8ipp566Qd8bgHI0ToHNbvvtt895552X73//+0mSkSNH5sc//nHzTeu8efOSvLFf6JupgjfTqO/kxRdfzH333Zd//OMfWbx4cRYvXpxTTjkl1157bUaMGJGrrroqy5YtS/LGnlAf/vCH8/DDD+fvf/9782dJ0r179+brv5koWJd3qm3kyJH5yU9+kuSNvVtffvnlJG/cPH/ta1/L4MGD1/kghdGjR+cXv/hFkuSGG27IAQccsNYepkma061J8ve//z1Lly7NjjvumPe///259dZb09TUlJdffjlz585Nz5491/vzAgDYEA888EAefvjhJG80NRcsWPCuq4HevK9qamrKH//4xw263siRI3P55Zc375H6wAMPNP/6TVOnTs1VV13VfO83ffr0XH311enatWu6deuWm2++OckbadWmpqZ13qftueeeWbhwYdasWZPHHnss8+fPX2c9y5Ytyw477JBtttkmd955Z/M95LruM5Nkn332yauvvpqpU6fmuOOOe9t8+++/f+bPn5+nnnoqy5cvb96rtdKpp56axsbGfO5zn2v+bPz48c33vFdddVU+/vGP5+KLL96gny8A5WicAm1i9OjRWb58ee6888586Utfak5K9unTJ5MnT06SnHHGGTnttNMycODAvOc971nvfNOmTcvhhx+e2tr//8faUUcdlalTp6Zv377513/91wwdOjT9+/fPRRddlI4dO+bKK6/M2LFj079//+YN/D//+c/nmmuuyYABA/L888+/4/XeqbZvfOMbefDBB9OvX78MGjQoDz74YJI3bp532WWXd1ymP3r06HTt2jU9e/bMN7/5zZx33nlJkjlz5uQLX/hCkjf2N+3fv38aGhqaUwy1tbX51Kc+lW7duqVfv34ZPHhwJkyY8LY9VAEAWmL58uU54YQT0qdPn/Tt2zdNTU057bTT1jtm0qRJGTt2bIYOHbrONOX6HHnkkTnssMMyaNCg9O3bN+PHj29OsiZvLNOfOXNmRo4c2fxZY2NjHn300TzzzDO58sor861vfSv9+/fP4YcfntWrV6/zPq179+756Ec/mj59+uSMM85I375911nPYYcdliVLlqRPnz654IIL0r9//yRZ533mm4477rgMHTp0nfex9fX1+f73v59DDjkkDQ0N+epXv9q82ujII4/Mk08+mQULFuSyyy7LjTfemIaGhjQ0NGTRokUb9HMEoGVqisrHNgOwyS1ZsiTDhg3LggUL1mruAgCwdTr++OPzhS98Ya3mLgBbFv/1DtDKbrjhhgwYMCCTJk3SNAUAqAJ9+vTJypUr89GPfrStSwGgBSROAQAAAAAqiD4BAAAAAO3Ok08+mYEDB6ZTp07ND0t+04IFC3LQQQflwAMPfMeH+7WUxCkAAAAA0O6sXLkyK1asyNFHH53p06envr6++djRRx+diy66KLW1tTn11FNz7bXXbvLr17/7KQAAAAAAm1enTp3SqVOndR578cUX84EPfCBJsnTp0la5ftU0Trt0rE+3Lh3bugzaqTWvC17zzpatXNXWJdCOdaiz6w3r9trq17P69aa2LmOr1rG+Np07Vs3tLBuoJjVtXQLtWI3fHsBGquu8XZ5//vm2LqNNHT64V55/6dUWzbGi9r3p3Llz8/tx48Zl3LhxGzRHU9P/v9durQX1VXOn2a1Lx5wybJ+2LoN26vnlK9u6BNqxvy78R1uXQDu2c9fO734SVWnBEy+0dQlbvc4d6zOk1y5tXQbtlL/YYn3qa/3+4J3V1ems884WN+3Q1iW0uedfejVzfr5hTc5Kjf/6x8yZM6dFc9S85W/Balvpz/WqaZwCAAAAAFuHHXbYIU888URqa2vTtWvXVrmGv2YDAAAAAMoripa9Slq9enVGjhyZe++9N4cddlhuueWWnHPOOUmSb3/72/nMZz6TT33qU5k8eXKrfE2JUwAAAACg3enQoUOmT5++1meHHHJIkmS//fbL7bff3qrXlzgFAAAAAKigcQoAAAAAUMFSfQAAAACgnCIbtE/plkziFAAAAACggsQpAAAAAFBedQROJU4BAAAAACpJnAIAAAAA5dnjFAAAAACgOmmcAgAAAABU0DgFAAAAAKhgj1MAAAAAoKTCHqcAAAAAANVK4xQAAAAAoIKl+gAAAABAedWxUl/iFAAAAACgksQpAAAAAFBOEQ+HAgAAAACoVhKnAAAAAEB51RE4lTgFAAAAAKgkcQoAAAAAbIDqiJxKnAIAAAAAVJA4BQAAAADKq47AqcQpAAAAAEAliVMAAAAAoLyiOiKnEqcAAAAAABUkTgEAAACA8qojcCpxCgAAAABQSeIUAAAAACiniD1OAQAAAACqlcYpAAAAAEAFjVMAAAAAgAoapwAAAAAAFTwcCgAAAAAoqfBwKAAAAACAaiVxCgAAAACUVx2BU4lTAAAAAIBKEqcAAAAAQHn2OAUAAAAAqE4apwAAAAAAFTROAQAAAAAq2OMUAAAAACjPHqcAAAAAANVJ4hQAAAAAKKf4v1cVkDgFAAAAAKggcQoAAAAAlGePUwAAAACA6qRxCgAAAABQwVJ9AAAAAKC86lipL3EKAAAAAFBJ4hQAAAAA2ADVETmVOAUAAAAAqCBxCgAAAACUVx2BU4lTAAAAAIBKEqcAAAAAQElFUlRH5FTiFAAAAACggsQpAAAAAFBOEXucAgAAAABUK4lTAAAAAGADVEfkVOIUAAAAAKCCxCkAAAAAUF51BE5bL3E6Y8aMTJo0aZPOuc8++2T48OH56Ec/mlNPPTUvv/zyJp0fAAAAACDZwpbq77zzzrn55pvz17/+NQcccEC+8Y1vtHVJAAAAAFBdiqJlry1EqzdOzz///Bx44IEZMWJEHn/88SRJv379MmbMmPTv3z/z5s1LkkyaNCkHH3xwTjvttJx88snvOu9JJ53UPBYAAAAAYFNq1cbp008/nZtuuim33357Jk+enHPPPTdJ8uyzz+ZnP/tZLr744kyZMiVPPfVU5s6dm1tvvTUHHXRQ6flramrWe/ySSy5JY2NjGhsb8+qqNS36LgAAtL233t+tWtPU1uUAALAVa9XG6eLFi7PffvslSRobG/PQQw8lSXr16pVOnTpljz32yNKlS/PYY4+lb9++SZKGhobS8xfvEu0dN25c5syZkzlz5qRLR8/BAgDY0r31/q5j/Ra16xQAwNajaOFrC9Gqd5s9evTIvffemySZM2dOevbsmWTtpGhRFOnevXsWLVqUJJk/f36puX/xi19k4MCBm7hiAAAAAICkVWOYu+66a3r27JmhQ4emY8eOmTJlyjrP22233dLQ0JBhw4ald+/e6dChwzrPe+655zJ8+PDU1tZmn332yfnnn9+a5QMAAAAAlbagBzy1RKs1Tg899NAceuihSZKzzjprrWO33XZbkjcSqVdccUWS5Oyzz059fX1+85vf5JFHHlnnnA888EBrlQsAAAAA0KzdbPw5ceLEzJo1K3V1dZk6dWrOOuuszJo1q/n4+PHj85nPfKYNKwQAAAAAqkW7aZxWLrs/99xz26gSAAAAAKDatZvGKQAAAADQzhWpmj1Oa9u6AAAAAACA9kbiFAAAAAAorzoCpxKnAAAAAACVJE4BAAAAgJIKe5wCAAAAAFQrjVMAAAAAgAoapwAAAAAAFexxCgAAAACUZ49TAAAAAIDqpHEKAAAAAFDBUn0AAAAAoLzqWKkvcQoAAAAAUEniFAAAAAAop4iHQwEAAAAAVCuJUwAAAACgvOoInEqcAgAAAABUkjgFAAAAADZAdUROJU4BAAAAACpInAIAAAAA5VVH4FTiFAAAAACgksQpAAAAAFBeUR2RU4lTAAAAAIAKEqcAAAAAQHnVETiVOAUAAAAAqCRxCgAAAACUUxT2OAUAAAAAqFYapwAAAAAAFTROAQAAAAAqaJwCAAAAAFTQOAUAAAAASiuKokWvDTFhwoQMGzYsp59++lqfX3311dl///1zwAEH5Nprr92UX6+ZxikAAAAA0O7MnTs3y5cvz8yZM7Nq1arcddddzccuuOCCzJgxIzNmzMgPf/jDVrm+xikAAAAAUFpRtOxV1uzZszNq1KgkyciRIzNr1qzmYz179swrr7yS5cuXp2vXrpv6KyZJ6ltlVgAAAACAFli6dGk++MEPJkm6deuWhQsXNh87+uijM2DAgBRFkZ///Oetcn2NUwAAAACgtA3dp7TSc889l8bGxub348aNy7hx4952Xrdu3bJs2bIkybJly7L99ts3H5s8eXIWLVqUJDnyyCPzsY99rEU1rYvGKQAAAACw2ey8886ZM2fOu543ZMiQ/PSnP82nP/3pTJ8+PSeffHLzsW222SZdunRJTU1NVq1a1Sp12uMUAAAAACil2ASvsgYOHJhOnTpl2LBhqaury5577plzzjknSTJ+/PgceOCBGTp06DrTqpuCxCkAAAAA0C5deOGFa72fOHFikuTkk09eK4HaGjROAQAAAIByipbvcbqlsFQfAAAAAKCCxCkAAAAAUFpTdQROJU4BAAAAACpJnAIAAAAApVXJFqcSpwAAAAAAlTROAQAAAAAqWKoPAAAAAJRSJCmqZK2+xCkAAAAAQAWJUwAAAACgtOrIm0qcAgAAAAC8jcQpAAAAAFBSkSZ7nAIAAAAAVCeJUwAAAACgnCKpksCpxCkAAAAAQCWJUwAAAACglCJJUSWRU4lTAAAAAIAKEqcAAAAAQGnVkTeVOAUAAAAAeBuJUwAAAACgtCZ7nAIAAAAAVCeJUwAAAACglCJJlQROJU4BAAAAACpVTeL0tdWv56FnlrV1GbRTnx3Sq61LoB2bMvPBti6Bdqxr545tXQLtVLXs+9TWampq2roE2im/N1ifjvUyRLyzV1etaesSaM/q2rqA9qGokntd/28BAAAAAFBB4xQAAAAAoELVLNUHAAAAAFqo8HAoAAAAAICqJXEKAAAAAJTWlOqInEqcAgAAAABUkDgFAAAAAEopYo9TAAAAAICqJXEKAAAAAJRUpKiSyKnEKQAAAABABYlTAAAAAKC0KgmcSpwCAAAAAFSSOAUAAAAASitSHZFTiVMAAAAAgAoSpwAAAABAKUWSpuoInEqcAgAAAABUkjgFAAAAAMopkqKojsipxCkAAAAAQAWNUwAAAACACpbqAwAAAAClVclKfYlTAAAAAIBKEqcAAAAAQGkeDgUAAAAAUKUkTgEAAACAUookTW1dxGYicQoAAAAAUEHiFAAAAAAozR6nAAAAAABVSuIUAAAAACitSgKnEqcAAAAAAJUkTgEAAACAcorCHqcAAAAAANVK4hQAAAAAKKVI0lQdgVOJUwAAAACAShKnAAAAAEBpRaojcipxCgAAAABQQeMUAAAAAKCCpfoAAAAAQGlFdazUlzgFAAAAAKgkcQoAAAAAlFZUSeRU4hQAAAAAoILEKQAAAABQSpGkqToCpxKnAAAAAACVJE4BAAAAgNKKVEfkVOIUAAAAAKCCxCkAAAAAUE6RFNUROJU4BQAAAACoJHEKAAAAAJRWVEnkVOIUAAAAAKCCxCkAAAAAUFqVBE4lTgEAAAAAKkmcAgAAAAClFCnSVCWRU4lTAAAAAIAKGqcAAAAAABUs1QcAAAAASquOhfoSpwAAAAAAbyNxCgAAAACUVng4FAAAAABAdZI4BQAAAADKKZIqCZxKnAIAAAAAVJI4BQAAAABKKZI0VUnkVOIUAAAAAKCCxCkAAAAAUFp15E0lTgEAAAAA3kbiFAAAAAAorbDHacvNmDEjkyZN2qRz7rPPPjn00ENz6KGH5tRTT92kcwMAAAAAJFtg4nTnnXfOjBkz2roMAAAAAKhKVRI43Tx7nJ5//vk58MADM2LEiDz++ONJkn79+mXMmDHp379/5s2blySZNGlSDj744Jx22mk5+eSTN0dpAAAAAABv0+qN06effjo33XRTbr/99kyePDnnnntukuTZZ5/Nz372s1x88cWZMmVKnnrqqcydOze33nprDjrooHec77nnnmteqj958uT1XvuSSy5JY2NjGhsb89qapk36vQAA2Pzeen+3yv0dAMBmV6RIU9Gy15ai1ZfqL168OAMGDEiSNDY25tvf/naSpFevXunUqVP22GOPLF26NI899lj69u2bJGloaMj//M//rHO+DVmqP27cuIwbNy5JssO227TwmwAA0Nbeen/XrUvHNq4GAICtWasnTnv06JF77703STJnzpz07NkzSVJTU9N8TlEU6d69exYtWpQkmT9/fmuXBQAAAADwjlo9cbrrrrumZ8+eGTp0aDp27JgpU6as87zddtstDQ0NGTZsWHr37p0OHTqs87w3l+onyR577JFf/vKXrVU6AAAAAPBWRfU8HKpVG6dv7kWaJGedddZax2677bYkbyRSr7jiiiTJ2Wefnfr6+vzmN7/JI488ss45H3jggVarFwAAAAAg2QyJ0w0xceLEzJo1K3V1dZk6dWrOOuuszJo1q/n4+PHj85nPfKYNKwQAAACA6lakOiKn7apxev7556/1/txzz22jSgAAAACAatauGqcAAAAAQPtWLXuc1rZ1AQAAAAAA7Y3EKQAAAABQSpGkqUoipxKnAAAAAEC7NGHChAwbNiynn376Wp+/8MIL+fSnP50RI0bknHPOaZVrS5wCAAAAAKVtrsDp3Llzs3z58sycOTPjx4/PXXfdlUGDBiVJvv3tb2fy5Mn58Ic/3GrXlzgFAAAAANqd2bNnZ9SoUUmSkSNHZtasWc3HFixYkO9+97sZPnz4Wp9vShKnAAAAAEBpRTZP5HTp0qX54Ac/mCTp1q1bFi5c2HzsjjvuyNy5c7PDDjvk2GOPzW233bbJr69xCgAAAABsNs8991waGxub348bNy7jxo1723ndunXLsmXLkiTLli3L9ttv33xs7733zr777pskqa1tnUX1GqcAAAAAQGkt3eN05513zpw5c971vCFDhuSnP/1pPv3pT2f69Ok5+eSTm4/tvffeeeqpp9K1a9esWbOmZQW9A3ucAgAAAADtzsCBA9OpU6cMGzYsdXV12XPPPXPOOeckeePhUJ/97GczYsSITJo0qVWuL3EKAAAAAJRSFEnR0sjpBrjwwgvXej9x4sQkSe/evTNjxoxWvbbEKQAAAABABY1TAAAAAIAKluoDAAAAAKU1bb6V+m1K4hQAAAAAoILEKQAAAABQUrFZHw7VliROAQAAAAAqSJwCAAAAAKVVR95U4hQAAAAA4G0kTgEAAACAUorEHqcAAAAAANVK4hQAAAAAKK2pOgKnEqcAAAAAAJUkTgEAAACAcgp7nAIAAAAAbDVWr169QedrnAIAAAAApRVFy16b05e+9KUkyZQpU9KvX79MmDCh9FiNUwAAAABgq/TQQw8lSf7yl7/k/vvvz9133116rD1OAQAAAIBSiiRFtpw9TpcvX55p06Zlp5122uCxEqcAAAAAwFbpP//zP3PPPfdk4sSJWbFiRY499tjSYyVOAQAAAICtygsvvJAk6dWrV04//fQkyWfsSc0AACAASURBVIoVK3LiiSeWnkPjFAAAAAAorWkLWKl/7LHHpqamJkXF06hqampy0003lZpD4xQAAAAA2KrcfPPNLZ7DHqcAAAAAQGlFUbTotTnNmTMnI0aMyP777581a9bkrLPOKj1W4xQAAAAA2Cp95StfydSpU7Ptttumvr4+f/vb30qPtVQfAAAAAChtM4dGW6S+vj477bRTampqkmSDEq8SpwAAAADAVqmxsTFnnHFGnn322UycODGDBw8uPVbiFAAAAAAopfi/f7YU3/ve93LjjTdm1113TZ8+fXLEEUeUHitxCgAAAABslYqiyMqVK5v/t6mpqfRYjVMAAAAAoJwiaWrha3P653/+59xyyy3ZY489MnPmzIwZM6b0WEv1AQAAAICt0jPPPJNf/epXSZIxY8ZkxIgRpcdqnAIAAAAApW3Ik+nbygsvvJAk6d27d6677roMGjQoc+bMyYABA0rPoXEKAAAAAGxVjj322NTU1KQoiixYsKD585qamtJzaJwCAAAAAKVtAYHT3HzzzS2eQ+MUAAAAANgqPfzww/nhD3+YJ554onmLgeuuu67U2NrWLAwAAAAA2HoUeWOP05a8NqeTTjopn/zkJ/Pss8/my1/+cj70oQ+VHqtxCgAAAABslTp27JiPfexj6dixYw477LAsWrSo9FiNUwAAAABgq9ShQ4esWLEiH/zgB/OVr3wlzz//fOmxGqcAAAAAQGlFC1+b05/+9Kd07tw5F198cYYNG1Z6f9PEw6EAAAAAgK3MNddcs87PZ82alWOOOabUHBqnAAAAAEBpTZv5AU8b47777lvn5zU1NRqnAAAAAEB1+ta3vtXiOTROAQAAAIByimQLCJxuEh4OBQAAAABsVV5++eUWz6FxCgAAAACUUiQpiqJFr83hqKOOSpKMGTNmo+dY71L90047LTU1Nes8dtFFF230RQEAAAAAWktRFBk/fnxmzpyZf/u3f1vrWNm+5nobp8cdd9zGVwcAAAAAbHW2hC1Ob7jhhtxzzz2ZPXt2jj322I2aY72N00MOOWSt9y+++OJmi9MCAAAAAGyMTp06ZciQIbnnnnuyZMmSPProo9lrr72y4447lp6j1B6nl112Wfr06ZOePXvmIx/5SEaNGrXRRQMAAAAAW6oiTUXLXpvTpZdemiOOOCI//elPc8QRR+SSSy4pPbZU4/QnP/lJ5s2bl/322y8PPPBA+vfvv9HFAgAAAABsDj/72c9yxx135NJLL83tt9+en//856XHrnep/pu6deuWDh06pCiK1NbWZuHChRtdLAAAAACwZSqSbEk7eRZFkdWrV6e+vj6rV69OU1NT6bGlGqejR4/OihUrctJJJ2XffffNYYcdttHFAgAAAABsDl//+tfT2NiY97///fnHP/6R73znO6XH1hQb8LSnV155Jdtuu+1GFdnWtuvUMQN77NTWZdBOvbxyVVuXQDs2pNcubV0C7djLK1e3dQm0U9fPezxLlq9s6zK2att32SaHfni3ti6DdqpjfV1bl0A79rw/n1mPnd7Tqa1LoB17pOm9mTNnTluX0aZ67LxdvvnJj7RojovnLd+sP8empqY8//zz2WmnnVJbW2rn0iQl9zi94YYbMmjQoOy///5Zs2ZNTjnllI0uFAAAAABgc6mtrc373ve+DWqaJiUbp+ecc05uu+22vO9970t9fX0ee+yxjSoSAAAAAGBLUGqP07q6umyzzTapqalJkrz++uutWhQAAAAA0A4VW9bDoVqiVOP0mGOOyac+9ak8/vjjOfHEE3Pccce1dl0AAAAAAC0yderUfP/7389rr72WoihSU1OT+fPnlxpbqnH67//+7znssMOycOHC7L333vnRj37UooIBAAAAgC1TU7acyOl3v/vd/PWvf82OO+64wWPXu8fpr3/963zoQx9Kjx49cs899+S6667L17/+dYlTAAAAAKDd69OnT7bffvuNGrvexOmFF16Y+fPn56WXXso+++yTadOmZfjw4Rt1IQAAAABgy1Zky9rj9JFHHsnee++d3r17J0lqampy3XXXlRq73sZply5d0rlz53Tu3DkDBgzQNAUAAAAAthhXXXXVRo9db+N0/vz5+cQnPpGiKLJgwYLmX29IZxYAAAAA2HoUW1DkdJdddsmll16aBx54IHvvvXe++MUvlh673sbp3Xff3eLiAAAAAADawtixY7P//vvns5/9bO64446ceOKJ+e1vf1tq7Hobp927d98kBQIAAAAAW4ctKHCa5557LmeccUaS5MADD9ygrUhrW6soAAAAAIC2tM0222TatGl56qmncs0116RDhw6lx2qcAgAAAAClNaVo0WtzmjJlSmbMmJHPf/7zufXWWzNlypTSY9e7VB8AAAAAYEu1yy675Ec/+tFGjdU4BQAAAABKKbJl7HH6+c9/Ppdffnn69euXmpqaJElRFKmpqcn8+fNLzaFxCgAAAABsVS6//PIkye23356uXbs2f758+fLSc9jjFAAAAAAopyhStPC1OR111FFrvf/c5z5XeqzEKQAAAACwVZk1a1buuOOO/O///m9++MMfJklWr16dJ598svQcGqcAAAAAwFalS5cu2WmnnbLddttlxx13TJJ07Ngxv/zlL0vPoXEKAAAAAJS2JTwcqn///unfv3/222+/9O/fP7W1tSmKIvPmzUv37t1LzWGPUwAAAABgq3TGGWektvaNFmhNTU3OPPPM0mMlTgEAAACA0jb3A55a4rXXXlvr/YoVK0qP1TgFAAAAALZKgwcPzhe+8IUMHz48M2bMyODBg0uP1TgFAAAAAEopkjS1dREb4Ac/+EGuv/76LFq0KEcffXSOPPLI0mM1TgEAAACArdagQYOy++67pyiKzJ07NwMHDiw1TuMUAAAAAChtS9rj9Oyzz8706dPzyCOPZK+99sq2226bP//5z6XG1rZybQAAAAAAbeLGG2/Mbbfdlr333ju33357dtppp9JjJU4BAAAAgNK2oMBpunTpkiSpra3N0qVLs2jRotJjJU4BAAAAgK3SF7/4xaxYsSITJkzI8OHDM2bMmNJjJU4BAAAAgHKKLWeP06ampjz22GPp3LlzPv7xj+fjH//4Bo2XOAUAAAAAtjq1tbW5++67N3q8xCkAAAAAUEqRpGnLCJwmSVasWJFBgwZl//33T11dXZLkoosuKjVW4xQAAAAA2CqdeeaZGz1W4xQAAAAAKK1I+4+cXnnllRk7dmwOOeSQLFy4MH369NngOexxCgAAAABsVa644ormX5922mkbNYfEKQAAAABQWtH+A6ebhMYpAAAAALBVmT9/fj7xiU+kKIq1fl1TU5Prrruu1BwapwAAAADAVuXuu+9u8RwapwAAAABASUWKLWCtfvfu3Vs8h4dDAQAAAABUkDgFAAAAAEopkjS1/8DpJiFxCgAAAABQQeIUAAAAACitSHVETiVOAQAAAAAqSJwCAAAAAOUUSVEdgVOJUwAAAACAShKnAAAAAEBpRZVETiVOAQAAAAAqSJwCAAAAAKUUSZqqI3AqcQoAAAAAUEnjFAAAAAAorSiKFr02xIQJEzJs2LCcfvrpbzu2YsWK7Lrrrpk+ffqm+mpr0TgFAAAAANqduXPnZvny5Zk5c2ZWrVqVu+66a63jl112Wfr169dq19c4BQAAAABKK1r4Kmv27NkZNWpUkmTkyJGZNWtW87FVq1Zl9uzZOfDAA1v6dd6RxikAAAAAsNk899xzaWxsbH5dcskl6zxv6dKl6dq1a5KkW7duWbp0afOxK664IieccEKr1lnfqrMDAAAAALzFzjvvnDlz5rzred26dcuyZcuSJMuWLcv222+fJFmzZk3+9Kc/5Xe/+13+9re/tVqdEqcAAAAAQDnF5ns41JAhQ/LXv/41STJ9+vQMHjw4SfLMM8/k8ccfz+GHH57//u//zllnnZUXX3xxk39VjVMAAAAAoN0ZOHBgOnXqlGHDhqWuri577rlnzjnnnOyxxx656667cuONN+aEE07Iueeem/e+972b/PqW6gMAAAAApW1AaLTFLrzwwrXeT5w4ca33Z599dqtdW+IUAAAAAKCCxCkAAAAAUEqRIk0tjpzWbJJaWpvEKQAAAABABYlTAAAAAKC0zbjFaZuSOAUAAAAAqNAqjdMZM2Zk0qRJm3TOvfbaKz/4wQ+a3/ft27dVn5oFAAAAALxdURQtem0ptpjE6e67757Zs2cnSR588MF06dKljSsCAAAAALZWrdo4Pf/883PggQdmxIgRefzxx5Mk/fr1y5gxY9K/f//MmzcvSTJp0qQcfPDBOe2003LyySevc66amprssssuefrpp3PNNdfk6KOPbs3SAQAAAIB1KIqWvbYUrdY4ffrpp3PTTTfl9ttvz+TJk3PuuecmSZ599tn87Gc/y8UXX5wpU6bkqaeeyty5c3PrrbfmoIMOWu+cn/zkJzNt2rT87W9/y+DBg1urdAAAAACgyrVa43Tx4sXZb7/9kiSNjY156KGHkiS9evVKp06dsscee2Tp0qV57LHH0rdv3yRJQ0PDeuccPnx4fv3rX2ePPfZITU3Nu9ZwySWXpLGxMY2NjVn9elMLvxEAAG3trfd3q9a83tblAABUnaJImoqiRa8tRX1rTdyjR4/ce++9SZI5c+akZ8+eSbJWw7MoinTv3j2LFi1KksyfP3+9c3bo0CGjR4/OIYcckldfffVdaxg3blzGjRuXJNmuU8eN+h4AALQfb72/277LNm1cDQAAW7NWa5zuuuuu6dmzZ4YOHZqOHTtmypQp6zxvt912S0NDQ4YNG5bevXunQ4cO6533zDPPTJLMmDFjU5cMAAAAALyLLSg02iI1RdH2X3XNmjWpr6/Pb37zmzzyyCM566yzNvk1tuvUMQN77LTJ52Xr8PLKVW1dAu3YkF67tHUJtGMvr1zd1iXQTl0/7/EsWb6yrcvYqm3fZZsc+uHd2roM2qmO9XVtXQLt2PP+fGY9dnpPp7YugXbskab3Zs6cOW1dRpvatVuXnHjg3i2a4+Zn67eIn2OrJU43xMSJEzNr1qzU1dVl6tSpOeusszJr1qzm4+PHj89nPvOZNqwQAAAAAKgm7aJxev7556/1/txzz22jSgAAAACA9SnS5gvYN4vati4AAAAAAKC9aReJUwAAAABgy9D2T0zaPCROAQAAAAAqSJwCAAAAAKUUSZqqJHIqcQoAAAAAUEHiFAAAAAAorUoCpxKnAAAAAACVJE4BAAAAgJKKFKmOyKnEKQAAAABABYlTAAAAAKCcwh6nAAAAAABVS+IUAAAAACilSNJUJZFTiVMAAAAAgAoSpwAAAABAaVUSOJU4BQAAAACopHEKAAAAAFDBUn0AAAAAoLSiStbqS5wCAAAAAFSQOAUAAAAASquOvKnEKQAAAADA20icAgAAAAClFLHHKQAAAABA1ZI4BQAAAADKKZKm6gicSpwCAAAAAFSSOAUAAAAASrPHKQAAAABAlZI4BQAAAABKq5LAqcQpAAAAAEAliVMAAAAAoJTi//6pBhKnAAAAAAAVJE4BAAAAgNKaqiNwKnEKAAAAAFBJ4xQAAAAAoIKl+gAAAABAaUVRHWv1JU4BAAAAACpInAIAAAAA5RRJlQROJU4B4P+1d/+xVtf1H8CflwtI9v0Khih0SzHYaAyR8DZAwO4lWVu1Of8oNmvF5oLhhuiqFQNcsvzB5kpWWdmmMXVO1qwt16yZXjW4lIwhMmbOKTicJPT1RpjKF+7n+8f36/3GUbgf5J57zr2fx2M7m4fP4X1eH3yfe9/3dZ/vzwcAAABqSZwCAAAAAKUUSYpUI3IqcQoAAAAAUEPiFAAAAAAorbcagVOJUwAAAACAWhKnAAAAAEBpRVGNyKnEKQAAAABADYlTAAAAAKC0igROJU4BAAAAAGpJnAIAAAAApRRxjVMAAAAAgMqSOAUAAAAASuttdAGDROIUAAAAAKCGxikAAAAAQA1b9QEAAACAcorCzaEAAAAAAKpK4hQAAAAAKK0igVOJUwAAAACAWhKnAAAAAEApReIapwAAAAAAVSVxCgAAAACUVo28qcQpAAAAAMB7SJwCAAAAAKX1usYpAAAAAEA1SZwCAAAAAKUUSc44cNoyEJXUn8QpAAAAAECNyiROR7QkY0a1NroMmtT0tgsaXQJN7MA/3mp0CTSx/f91pNEl0KTe/u/jjS4BKu3oMZ9BTu7s0ZX5UZgP4Fhvb6NLgOZWJMWZRk4lTgEAAAAAhia/ZgMAAAAASjvja5wOERKnAAAAAAA1NE4BAAAAAGrYqg8AAAAAlNabauzVlzgFAAAAAKghcQoAAAAAlFLEzaEAAAAAACpL4hQAAAAAKKlIUZHIqcQpAAAAAEANiVMAAAAAoLSKBE4lTgEAAAAAakmcAgAAAACl9aYakVOJUwAAAACAGhKnAAAAAEApReEapwAAAAAAlSVxCgAAAACUVlQkcipxCgAAAABQQ+IUAAAAACitIoFTiVMAAAAAgFoapwAAAAAANWzVBwAAAABK663IXn2JUwAAAACAGhqnAAAAAEApxQA8TseNN96YhQsXZtWqVSf8+fLlyzN//vwsWLAgu3bt+uAndAoapwAAAABA09mxY0eOHDmSp59+OkePHs0zzzzTd+y73/1utmzZknvvvTc333xzXd5f4xQAAAAAKK0oijN6lLVt27YsXrw4SXLllVemu7u779jFF1+cJBk1alRaW1sH9gT/j8YpAAAAADBoDh48mPb29r7H3Xff/b6v6+npyTnnnJMkGTt2bHp6et7zmtWrV+f666+vS50j6zIqAAAAADD8FMlphEbf14QJE7J9+/Z+Xzd27NgcPnw4SXL48OGMGzfuhON33nlnpk+fngULFpxZQSchcQoAAAAANJ158+blj3/8Y5Lksccey9y5c/uO/eEPf8jWrVuzdu3aur2/xikAAAAAUNKZXd/0dK5xOnv27IwZMyYLFy5Ma2trLrzwwtxyyy1JkpUrV+bll19OZ2dnli9fXpcztVUfAAAAAGhKGzduPOH5mjVrkiR//etf6/7eGqcAAAAAQClFkt4zvMbpUGGrPgAAAABADYlTAAAAAKC0ItWInEqcAgAAAADUkDgFAAAAAEorqhE4lTgFAAAAAKilcQoAAAAAUMNWfQAAAACgtKIie/UlTgEAAAAAakicAgAAAAClFEXSW43AqcQpAAAAAEAtiVMAAAAAoLQi1YicSpwCAAAAANSQOAUAAAAASiuqETiVOAUAAAAAqCVxCgAAAACUVlQkcipxCgAAAABQQ+IUAAAAACilSNJ7hoHToZLkHCp1AgAAAAAMGolTAAAAAKCkwjVOAQAAAACqSuIUAAAAACitGnlTiVMAAAAAgPfQOAUAAAAAqGGrPgAAAABQThE3hwIAAAAAqCqJUwAAAACglCJJbzUCpxKnAAAAAAC1JE4BAAAAgNJc4xQAAAAAoKIkTgEAAACA0qqRN5U4BQAAAAB4j7o2Tru6urJ27doBHXPatGnp7OzMnDlzsnXr1gEdGwAAAAA4taIozugxVAy5rfoTJkzIE088kVdffTXXXnttHn300UaXBAAAAAAMM4OyVX/Dhg2ZP39+Fi1alFdeeSVJcskll+Saa67JpZdemp07dyZJ1q5dmyuuuCIrV67M0qVLTzlmW1tb3nnnnXqXDgAAAAD8m6I4s8dQUffG6YEDB/L4449ny5YtWb9+fW677bYkyeuvv5577rknd911VzZt2pTXXnstO3bsyFNPPZUFCxb0O+7zzz+fcePG1bt8AAAAAKCC6t443bt3b2bOnJkkaW9vz4svvpgkmTp1asaMGZO2trb09PRk3759mTFjRpJk1qxZJx3v4MGD6ezszA033JBbb731lO999913p729Pe3t7Tl6vHeAzggAgEY5YX137HijywEAqJyiSHqL4oweQ0Xdr3E6efLkPPvss0mS7du3Z8qUKUmSlpaWvtcURZGLLrooe/bsSZLs2rXrpOO9e43TMpYtW5Zly5YlScZ+aPQHqh8AgObx7+u7cWef1eBqAAAYzureOJ04cWKmTJmSyy+/PKNHj86mTZve93WTJk3KrFmzsnDhwkyfPj2jRo2qd2kAAAAAwGkaQqHRM9JSFM1zqseOHcvIkSPz0EMP5aWXXsrq1asHbOyxHxqduVMvGLDxGF4+cf45jS6BJvb64bcaXQJNbP9/HWl0CTSp3fvfyJvv/HejyxjWxp19Vjo+OanRZQBD0LHepvkxmCY0euSg3EebIeqVjM/27dsbXUZDjRk1MheO/48zGuOcj00dEv+OdU+cno41a9aku7s7ra2t2bx5c1avXp3u7u6+4ytWrMiSJUsaWCEAAAAAVFmRItX4BVRTNU43bNhwwvPbbrutQZUAAAAAAFUmfw4AAAAAUKOpEqcAAAAAQHNrnjsm1ZfEKQAAAABADYlTAAAAAKCUIklvRSKnEqcAAAAAADUkTgEAAACA0ioSOJU4BQAAAACoJXEKAAAAAJRWpBqRU4lTAAAAAIAaEqcAAAAAQDmFa5wCAAAAAFSWxCkAAAAAUEqRpLcikVOJUwAAAACAGhKnAAAAAEBpFQmcSpwCAAAAANSSOAUAAAAASisqEjmVOAUAAAAAqKFxCgAAAABQw1Z9AAAAAKC0amzUlzgFAAAAAHgPiVMAAAAAoJQiRXrdHAoAAAAAoJokTgEAAACAcoqkIoFTiVMAAAAAgFoSpwAAAABAaUVFIqcSpwAAAAAANSROAQAAAIDSKhI4lTgFAAAAAKglcQoAAAAAlFIkKVKNyKnEKQAAAABADYlTAAAAAKC03moETiVOAQAAAABqSZwCAAAAAKUVRTUipxKnAAAAAAA1NE4BAAAAAGrYqg8AAAAAlFMkFdmpL3EKAAAAAFBL4hQAAAAAKKVIUqQakVOJUwAAAACAGhKnAAAAAEBpvdUInEqcAgAAAADUkjgFAAAAAEoqUhTViJxKnAIAAAAA1JA4BQAAAABKq0jgVOIUAAAAAKCWxCkAAAAAUEqRuMYpAAAAAEBVSZwCAAAAAKX1NrqAQSJxCgAAAABQQ+IUAAAAACincI1TAAAAAIDK0jgFAAAAAKhhqz4AAAAAUFpFdupLnAIAAAAA1JI4BQAAAABKKeLmUAAAAAAAlSVxCgAAAACU1tvoAgaJxCkAAAAAQA2NUwAAAACgtKIozuhxOm688cYsXLgwq1atOuHPd+/enQULFmT+/PnZtWvXQJ5eH41TAAAAAKDp7NixI0eOHMnTTz+do0eP5plnnuk7tm7dujz44IPZvHlz1q1bV5f31zgFAAAAAEorijN7lLVt27YsXrw4SXLllVemu7u779gbb7yRj3/842lra0tPT89An2ISjVMAAAAAoAn19PTknHPOSZKMHTv2hAZpb+//36LqdLf/lzWyLqM2oVEfPid/P+uCRpfRNA4ePJgJEyY0uoym8fd/NLqC5mFuvJ8PNbqApmF+vI9zP9LoCpqG+XGiD/2jKvcabZyRZ/9n9o84r9FlNA2fQU7F/KghQnQC84NTMT9OtHfv3kaX0HCf+9zncujQoTMa46233kp7e3vf82XLlmXZsmXved3YsWNz+PDhJMnhw4czbty4vmMtLS19/z1iRH2+sFemcXqm/0OHm/b29mzfvr3RZdCEzA1OxfzgVMwPBpv13Yl8BjkV84NTMT84FfODWo8++uigvde8efPy85//PF/+8pfz2GOPZenSpX3HPvKRj2T//v0ZMWJEXyp1oPk9GwAAAADQdGbPnp0xY8Zk4cKFaW1tzYUXXphbbrklSXLzzTdnyZIl+dKXvpT169fX5f0rkzgFAAAAAIaWjRs3nvB8zZo1SZKZM2dmy5YtdX3v1u9973vfq+s70LQuu+yyRpdAkzI3OBXzg1MxP6CxfAY5FfODUzE/OBXzg6pqKep12ykAAAAAgCHKNU4BAAAAAGponAIAAAAA1NA4BQAAAACooXE6xHV1dWXt2rUDOua0adPS2dmZz372s7nuuuvyz3/+c0DHZ3DVa450dHSko6Mj11133YCOzeCpx9y4+OKLc8cdd/Q9nzFjRtyDcOiq5/eYOXPmZOvWrQM6Ngwn1nj0xxqPk7HGoz/WeFDeyEYXQPOZMGFCnnjiiSTJpk2bsm7dutx5550NropmMmHChHR1dTW6DJrQRz/60Wzbti1J8sILL+Tss89ucEU0m3e/x7z66qu59tpr8+ijjza6JKgMazz6Y43HyVjj0R9rPIYridNhYsOGDZk/f34WLVqUV155JUlyySWX5Jprrsmll16anTt3JknWrl2bK664IitXrszSpUv7HffrX/96399laKvXHGHoG8i50dLSkgsuuCAHDhzIww8/nKuvvnqwToM6qsfXj7a2trzzzjv1Lh2GPGs8+mONx8lY49Efazzon8bpMHDgwIE8/vjj2bJlS9avX5/bbrstSfL666/nnnvuyV133ZVNmzbltddey44dO/LUU09lwYIFpcdvaWmpV+kMkoGeIwcPHuzbxrV+/frBOg3qoB5fP6666qr85je/yZ///OfMnTt3ME6DOqrX95jnn38+48aNq3f5MKRZ49EfazxOxhqP/ljjQTkap8PA3r17M3PmzCRJe3t7XnzxxSTJ1KlTM2bMmLS1taWnpyf79u3LjBkzkiSzZs0qPX5RFANfNINqoOfIu9u4urq6ctNNN9X/BKibenz96OzszIMPPpi2tjY/lA8DAz1HDh48mM7Oztxwww259dZb638CMIRZ49EfazxOxhqP/ljjQTmucToMTJ48Oc8++2ySZPv27ZkyZUqSE1MERVHkoosuyp49e5Iku3btKjX2fffdl9mzZw9wxQy2es4RhrZ6zI1Ro0bli1/8Yj7zmc/kX//6V50qZ7AM9Bz592ssAqdmjUd/rPE4GWs8+mONB+VonA4DEydOz+q2qQAABHdJREFUzJQpU3L55Zdn9OjR2bRp0/u+btKkSZk1a1YWLlyY6dOnZ9SoUe/7und/UzRixIhMmzYtGzZsqGf5DIJ6zJGOjo4k/3sNmwceeKBepVNnAz033vXtb387SdxgYhio1xwB+meNR3+s8TgZazz6Y40H5bQU9uhUyrFjxzJy5Mg89NBDeemll7J69epGl0STMUc4GXOD/pgj0Dg+f/THHOFkzA36Y45QZRqnFfOd73wn3d3daW1tzebNm/ODH/wg3d3dfcdXrFiRJUuWNLBCGs0c4WTMDfpjjkDj+PzRH3OEkzE36I85QpVpnAIAAAAA1BjR6AIAAAAAAJqNxikAAAAAQA2NUwAAAACAGhqnwLCyd+/ejB8/Ph0dHfn0pz+d3//+96c9RkdHRw4dOpSdO3dm48aNJ33dXXfddVrjzpgx47RrAQDAGg+AxtA4BYadefPmpaurK7/+9a/zzW9+s+/Pjx8/flrjzJo1K6tWrTrp8dNdVAMA8MFZ4wEw2EY2ugCAevnYxz6WN998M1OnTk1nZ2fefPPN/PKXv8zy5cuzb9++JMndd9+dqVOnZuPGjbn//vvziU98Ij09PUmSrq6u/OpXv8qPf/zjPPbYY7npppty1llnZfHixfnwhz+cl19+OR0dHVmxYkXmzp2bFStW5O2338748eNz3333ZfTo0fna176W/fv357LLLmvkPwUAwLBhjQfAYNE4BYat5557Luedd152796d9evXZ9KkSfnpT3+aT33qU7n33nvz3HPPZd26ddm4cWMeeOCBdHd358iRI7n44otPGKe3tzerVq3KU089lfHjx+f48eNpbW3NL37xi3R1dSVJlixZkttvvz0zZ87Mj370o2zatCnnn39+zj333Nx///158sknP9CWMgAATmSNB8Bg0TgFhp3u7u50dnZm1KhR+dnPfpavfvWrmTRpUpJk9+7d2bZtWx5++OEk6UsVXHrppRk5cmTGjRuXadOmnTDeoUOHMnHixIwfPz5J0tra+p733LNnT66//vokydtvv50vfOEL6enpydy5c5Mkc+bMqdv5AgBUgTUeAINN4xQYdubNm5dHHnmk7/m/L4KnT5+emTNnZvny5UmSo0ePpqenJ7t27crx48dz5MiRvPDCCyeMd9555+Vvf/tb3njjjZx77rnp7e3NiBEj0tLS0veaT37yk/n+97/ftyA/evRoHnnkkTz55JP5yle+kr/85S/1PGUAgGHPGg+AwebmUEClfOMb38gzzzyTRYsWpbOzMz/5yU9y/vnnZ8mSJZkzZ05WrFiRyZMnn/B3RowYkR/+8If5/Oc/n87Oztx+++1J/vfGAldffXV++9vf5o477si3vvWtLFq0KIsWLcqf/vSnXHXVVTl48GA6Ojryu9/9rgFnCwBQDdZ4ANRDS1EURaOLAAAAAABoJhKnAAAAAAA1NE4BAAAAAGponAIAAAAA1NA4BQAAAACooXEKAAAAAFBD4xQAAAAAoIbGKQAAAABADY1TAAAAAIAa/wP350CBsDz9nwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1656x720 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot confusion matrix\n",
    "plot_cm(np.concatenate(out['y_true']),np.concatenate(out['y_pred']),\n",
    "        np.concatenate(out['sy_true']),np.concatenate(out['sy_pred']),decoder+'_confusion')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f22a19a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAHiCAYAAABMVsECAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhcVZ3/8fcnIUACISAgIAmrGnAZZZeAIos4E0UIiLgBruMomzP6IPBzZ0RgHEU2d0RRVEDEBdzYBUGWAIPKJhDoIAFZBJKQhCTf3x/nFimKrurq7rr3VFd9Xs/Tz+2qe+rebyfpT+52zlFEYGaW07jcBZiZOYjMLDsHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXXN0EkabqkwySdKelWSUslhaRPjnK7e0i6SNIjkp6WdLukL0havVO1m/W6lXIXUKEPA0d0coOS/hP4MhDAH4CHgNcCxwD7Sdo5Ih7p5D7NelHfHBEBfwa+BLwL2BI4azQbk7QV8L/AMuBNEbFLRLwN2By4BJgOfH1UFZv1ib45IoqIb9e/lrR8lJs8GhDw3Yj4dd1+Fkp6P3AP6ahoi4i4fZT7Mutp/XRE1DGSVgbeVLw8u3F9RNwHXF28nFVVXWZjlYNoZF4KTCq+v6FJm9r7W5VfjtnY1jenZh22abH8Z0Q81aTNQEPbUZE0dYgm44F1gXnAvIhY2on9mlXBQTQyk4vlghZt5hfLNTq0z4GhmzxrGjC3Q/s1K52DqAcNDAwMJ7TMRkqd2pCDaGRqp2OrtWhTe6DxyQ7tc9oQ69cHru/Qvswq5SAamTnFck1Jk5tcJ5rW0HZUIqLlqZbUsf+czCrnu2YjcwewsPh+2yZtau/PLr8cs7HNQTQCEbEEuLB4+c7G9ZI2BmYUL39WVV1mY5WDqAVJhxadWL8/yOrjSX3M3ivpX+s+Mwn4Dul2+k/9VLXZ0PrmGpGkrYHT697avFh+SNKb696fFREPFt+vQ+ozNq9xexExW9LHSJ1eL5J0BfAwqdPrBqTTt//o7E9h1pv6JohIz/PsMMj7U4uvmlXa3WBEfEXSrcDHgO1Jd9HuB74IfLHFw45mVkcRkbsG64DiyesBgIGBAaZOHepBbLNR69itWl8jMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXnIDKz7BxEZpadg8jMsnMQmVl2DiIzy85BZGbZOYjMLDsHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXnIDKz7BxEZpadg8jMsnMQmVl2DiIzy85BZGbZOYjMLDsHkZll5yAys+wcRGaWnYPIzLLruyCStL+kyyU9LmmBpFskHSlpwgi2tZqkoyXdIOlJSc9ImifpV5LeUkb9Zr1opdwFVEnSScARwFLgUmA+sBtwArCXpD0j4uk2t7U2cCXwsmI7fwT+CbwYeBPwJkknR8QRHf9BzHpM3xwRSdqHFELzgR0i4o0RsR/wEuBWYGfg2GFs8tOkELoR2LjY3gERsQ0piJYCh0t6TSd/DrNe1DdBBBxTLI+PiNm1NyPiEeAjxctDJU1pc3u7FcsTIuKx+hURcRFwWfFyxxHWa9Y3+iKIJG0IbFe8PLtxfURcBQwAqwAz29zsojbbPdJmO7O+1RdBBGxVLB+LiHubtLmhoe1Qfl0sPyHpBfUrJM0EdgXmAb8YTqFm/ahfLlZvWizvb9FmoKHtUE4AtgfeCNwn6WpWXKzeBrgaeH9EPDH8cp9P0tQhmqzfif2Y5dAvQTS5WC5o0WZ+sVyjnQ1GxAJJewHHAR8jBVLNo8DFwAPDrLOVgaGbmI1N/XJq1nGSNiAd9RwGfBLYDFiddJR0I/AZ4CpJk5tuxMyA/jkieqpYrtaizerF8sk2t/k90gXwIyPif+rev17Sm0lh9Crg46RQGq1pQ6xfH7i+A/sxq1y/BNGcYtnql7m2bk6LNsCzd+HeULz8UeP6iHhG0nnAK4E96EAQRcTcIWoa7S7MsumXU7ObiuXakppdjN62WM5usr7eRnXfNzuCql2kfkGT9WZW6IsgKo4maqct72xcL2ln0hHRYuCiNjZZfxF6hyZtak9UN3tcwMwKfRFEheOK5VGStq69WfQZO714eWr97XZJsyTdLumS+g1FxP2sCLavStqkfr2kdwMHFC+f9wClmT1Xv1wjIiIukHQycDhwbREuC4DdgTVJd8A+1fCxKcB0YNVBNvk+UjeOLYHbJF1Leop6S+DlRZsfAD/s8I9i1nP6JogAIuKI4sHDQ4AZwATgbuB44CsRsWQY2/qzpFcA/wn8G+kO2irA48BvgTMi4pwO/whmPUkRkbsG64DiyesBgIGBAaZOHepBbLNR69it2n66RmRmXcpBZGbZOYjMLDsHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXnIDKz7BxEZpadg8jMsnMQmVl2DiIzy85BZGbZOYjMLDsHkZllV3oQSdpH0viy92NmY1cVR0TnA/dL+qKkl1SwPzMbY6oIovnABsCRwO2SLpd0oKRVK9i3mY0BVQTR+sD7gWsAAa8DzgQelHS6pG0qqMHMuljpQRQRCyPiuxGxM7AF8CXgYWAK8CHgOkk3SzpE0ppl12Nm3afSu2YRcWdEHAlMBWYBFwHLgX8BTgb+LukHknarsi4zyyvL7fuIWBYRP4+IvYCNgGOAu4FVgXcAv5f0N0mfkPSCHDWaWXW64Tmi1UmnaasDQbqOJGAz4DhgjqSP5yvPzMq2Uo6dSpoIvA14H7Bz7W1gDvAd4KfA7sB/AC8HTpC0JCJOrr5aMyubIqK6nUk7kMLnAGAyKXyWAr8Evgn8LhoKkvQJ4IvAnRGxRWXFjjGSpgIDAAMDA0ydOjVzRdYH1LENlR1EktYBDiTdwt+y9jZwD/Bt4LsR8dAQ23gcmBQRq5RZ61jmILIMOhZEVZyazQUmkIp+BrgA+FZEXDyMbTwBrFFCbWbWBaoIopWBu0hHP2dGxD9GsI0DSHfUzKwHVRFEu0XE5aPZQET8qUO1mFkXquLJ6svL3oeZjW2lHxFJmgzsCjwVEZcN0XY30vNEl0bE/LJrM7PuUMUDje8Gfgb8axtt9yvavqOsYiTtX4wA8LikBZJukXSkpAmj2Obekn4haZ6kJZIelvRHSZ/uZO1mvaqKINq3WP6wjbZnkO6u7VdGIZJOAs4BdgKuA35D6mJyAnBp8aDlcLa3sqRzSHcC9wD+ApwH/BnYHDi8c9Wb9a4qLla/FFgC3NpG29lF25d2ughJ+wBHkMZH2iUiZhfvrwNcSnrC+1hgON1JvgXsTwqiD0bEI3X7Gwds35nqzXpbFUdE6wHzG5+YHkzR5inSGEaddkyxPL4WQsU+HwE+Urw8VNKUdjYmaXfgINLRz9vqQ6jY7vKIuHb0ZZv1viqC6ClgTUlDPhVdjNq4JvB0JwuQtCGwXfHy7Mb1EXEV6ankVYCZbW72sGJ5UkQ8M+oizfpYFUH052I/b2mj7VuA8cBtHa5hq2L5WETc26TNDQ1tmyomA9i9eHmlpPUlfVTS1ySdJOlgSauPsmazvlHFNaLzgV2AL0v6U0TcP1gjSRsD/0saCuT8DtewabEcdN+FgYa2rWxGeswA4DXA6XWva/5H0tsj4tK2q2yh6EvWShmns2aVqOKI6JvA34ANgZslfVrStpLWLb62lfQZ4Kaizb2kX+xOmlwsF7RoU3tuqZ0+bWvXff8d4EbSqd9k4NWkkSfXBX7ewZlLBob4ur5D+zGrXOlHRBGxWNKbSbfKNwE+U3w1qo1H9KaIWFR2XaNU3+v4AeCNEbG4eH2LpLcANwOvAI4ijTxgZk1UMkJjRNxJuvZyAjCPFaMw1r7mkcYc2ioi7iihhKeK5Wot2tROrZ4cxvYgdeRdXL8yIpYB3yhe7tFWhUObNsTXds0/atbdKhuhMSKeAI4Gji6uB61XrJrX7LpRB80pltNatKmtm9OiTf32asPa3tOkTe39DdrY3pAiYm6r9VLHhoYxq1yWoWIj4j7gvgp3eVOxXFvSpk3unG1bLGcPsu45ImK+pDtI0yOt06RZ7X33mTMbQjcMnl+64miidjH3nY3rJe1MOiJaTLrQ3I5zi2WzU683FMvr2tyeWd+qPIiUvEDSNEkbNfsqYdfHFcujJG1dV8/arLhLd2pxCllbN0vS7ZIuGWR7JwOPAzMlfah+haS3A++qa2dmLVR2aiZpP1JXitcw9GiLQYdri4gLJJ1M6oh6bREuC0gPJq4JXA18quFjU4Dpg9UbEY9IOgD4BfB1SYeRHsTcnBUPRR4bEe0eYZn1rUqCSNLppOml272iWsqV14g4QtLVwCHADNJY2ncDxwNfiYglw9ze7yW9itSPbQ9gb9Jdt4uAr0bE7zpZv1mvqmIWj1mkecqWAp8FLiRdPP4HsCPp7tmupCOVcaRnbv6vuKBtbfIsHpZBxw4YqrhG9AHSqdanIuK4iLileH9ZRNwTEddExHHAv5DC6Tuk0DKzPlFFEG1TLL/dat8R8TDwYdJt72Mws75RRRCtRRqP6NG695YyyFPOEXElsBD4twrqMrMuUUUQPQosb3jvcWCSpBcM0l506GlkMxsbqgiiAWCNYjaPmr8Wy+cMqC9pBjCR1r3kzazHVBFEtQHHXlP33i9IRz6nSTpI0ssk7UsaYD+AKyqoy8y6RBVBdD4pdN5V997XSNNQTwG+SxpY/1xgY9IwsZ+toC4z6xJVBNGlpFEP/1/tjWK8oV1I4bOEFc8jXEOaorqdGT/MrEdUMTBaMEhP+4iYBxxQTGy4DmkmWPdUN+tDVUw5XRs0/4+NU+4AFDNgPFh2HWbWvaroa3YB6bmhwW7Vm5lVEkSPQRpMrIJ9mdkYVMXF6r8AUyS1MzuGmfWhqqYTGs+KmVHNzJ6jirtmP5S0PfC5Ykrpr0TEY2Xv18zGjirumtVmOl1I6lX/CUl/Iw35sazJxyIidm+yzsx6TBUXq18/yD63KL6aKXe0NjPrKlUE0ecq2IeZjWFVXCNyEJlZS30xr5mZdTcHkZll5yAys+yquH3f7BZ9KxERlU3+aGZ5VfHLXspkiWbWO6oIol2HWD8F2I40/9kEUleQv5ddlJl1j9Jnem2XpLWBi4HVgW0j4onMJY0pnunVMhhTM722pZj37MPA5sAnM5djZhXqmiACiIhrSVMJ7Zu7FjOrTlfdmZI0nlTTi3LXYp3z1re+lYULFwIwadIkzjvvvMwVWbfpqiACZgKrAPNyF2Kds3DhwmeDyGww2YNI0srAVGBv0rWhAH6XtSgzq1S3PdAo4GHgMyWVY2ZdqIqL1WrzazHwY2CHiHjePGhm1ru64YHGpcDjwF3FHGdm1meqGI/oirL3YWZjW1c9R2Rm/amSIJK0hqTV22i3uuc/M+s/pQeRpP1I14C+0Ubzs4DHJe1dblVm1k2qOCI6oFh+p4223yTdQXt7eeWYWbepIoheTZq/7Oo22l5atN2q1IrMrKtUEUQbAk9ExOKhGhZt/on7mpn1lSqC6BnSGENDkqR225pZ76giiO4BVpY0o422M0idXv1ktVkfqSKIfku6AH2ipAnNGhXrTiR1ev1NWcVI2l/S5ZIel7RA0i2SjmxV2zC2PVNSFF8Xd6Jes35QRRCdBMwHdgSulPS6+pVKdgH+ULRZCHyljEIknQScA+wEXEcKvI2AE4BLJU0cxbbXAr5FClIzG4bSgygiHgLeSbpWtD1wmaQnJP1F0l9IF6cvLdY9A7w7Ijo+eL6kfYAjSKG4Q0S8MSL2A14C3ArsDBw7il2cAqwHfH20tZr1m8oGz5e0Peno6DVNmlwDfDQiri9p/9eRZgv5ZER8oWHdzqQjssXAesMduF/SLOB84H+AvwLfBS6JiD06UXubNbQcPL9+lMSqNe530qRJWerw6JAd17HB8ysbGC0irgNmSHop6RRsvWLVPOCaiLirrH1L2pAUQgBnD1LbVZIGgGmkUSJ/NIxtr0M6CroD+DRd+jBmN42S2C11WPeofITGiLgTuLPi3dYekHwsIu5t0uYGUhBtxTCCCPgasA6wb0QsSk8gmNlwZB8qtiKbFsv7W7QZaGg7JElvB94KfDUi2nlyfMSKU69W1h/O9mLcqG8Stm95wzBTFe5bjfu2rlTFULGbkmZvvS8ivjpE24+RnsQ+KSJahcZwTS6WC1q0mV8s2+r9L2l94DTgbuCYkZfWtoGhm7Qnxk1g/jYHdmpzXW31G89yGI0BVdy+P4h0t6qdq+KrFm3Hwm/JN4G1gA9EhC96mI1CFUE0s1he0EbbH5CuxL+pwzU8VSxXa9Gm1rXkyaE2JulgYC/g6xFx+ehKa9u0Ib62a/5Rs+5WxTWijYEF7ZxqRcR9khYUn+mkOcVyWos2tXVzWrSpmVUst5N0ecO62rWaberWvT0iRjVXW0TMbbXeF8ltLKsiiNak9bWZRkuAtTtcw03Fcm1Jmza5c7ZtsZw9jO1u22LdmsAuxferDmObZn2nilOzR4A1JQ0ZLkWbNUkjOnZMcTRRe1DynYPsd2fSEdFi4KI2trdPRGiwL+C9RbNL6t6f05mfxKw3VRFE1xXL97XR9v2ka0Q3lFDHccXyKElb194swu/04uWp9U9VS5ol6XZJl5RQj5kVqgii75PC5fOtxqIu+oJ9nnR37fudLiIiLgBOJl2UvlbSryWdB/wNeCVpBMlPNXxsCjAd2LzT9ZjZClXMa3aBpAtJd8LOl3QlcCErxhzauFj3OlJg/SYizi2pliMkXQ0cQhr7aALpOaDjga9ExJIy9mtmrVXS6VXSaqRb87Ujosad1m75/Bw4MCLmY8MyVKfXXXfdldrfdUClTzdntfyZZ/9xSeKyyy7LWk676jspd3Fn3bHV6TUiFgCzJO0JvIdBOr0CZ0bE76uopx/V/4cjeH63iz5Q1UgTndBNnZSrUGlfs4j4HfC7odpJeoNDyax/dE2nV0kvIR0tHUiaxaNrausFknxq5oc+u1bWX3ZJk0kTML6HdLoG6cxh7BxDjxETJ05ccajfZ51ea6ehEyeOeCRgK1mWIJK0Byl89gEmsuKi153AecBPc9RlZnlUFkSSNmfFqdc0VoRPkB42/FFE/KWqevqZlj+TjhSq4vGIbAilBpGk1YG3kbo91OY1E/AP4MekcYoAToyIIXu9W+dk/QV1OFiDUoJI0m6ko59ZwCRS+CwCfgGcRXpocZmkw5puxMz6RseDSNK9pLnCahedryCFz3kR8VSrz1p5cs2cAd01i4d1pzKOiDYmBdA3gP+OiAdK2IcNU84nc2fOnPmcp4QvumjIAQ6sz5TV6VXAB4GzJX1Q0pol7cfMekAZR0Q7kob8eBvwWtIMqqdIuojU3+xX7lxari222IJx44b3f8wpp5zCwQcf/Lz3v/e973HYYSO7lPfkk4Pff9h99925/vrhz6N50EEHceqppz7v/dtuu40ddthh2NsDuPjii9l+++2f9/7RRx/NaaedNuztbbHFFlx33XWDrttwww156qn2rk7Un85Onz590DY/+9nPBv07a8fcuXNZY43nzxMxa9YsLrmkvVFvmv39jkTHgygi/gT8SdIRwL6ki9a7kZ4Z2ht4ohh+44ed3rclCxYMZ0DM5JlnBr+T9cwzz7T9y9OuhQsXjmibixYtGvT95cuXj7jGZcuWDfr+4sWLR7TNVn/2Tz311Ii2uXz58kHfX7p0adf83YxWaeMRRcSiiDg7IvYENiGN9XM3aQTGD5Dmu695RVl1mFn3q6r3/VzgC8AXJO1Eeq5of9J8YwL+UNxtO490d62MERr7xmqrrTbsU7MJEwZ/yHDChAlMnjx50HUjNWnSpBFtc9VVBx/6e9y4cSOucfz48YO+v8oqq4xom6ut1nyimOFsr/7UrNnf5UorrdQ1fzejVcl4RIPuWJpImiX1PcDreW4fs/siYrMshY1RQ41HlJPvmg3fGPkz61gv4iqGih1URDwdEWdFxO7AZsDnSFP5iM5PJ2RmXSxbENWLiPsi4nMRsTnpwnbHx6w2s+7VdWP+FDOnXp65DDOrUFccEZlZf3MQmVl2DiIzy85BZGbZdd3FarNuUT+3WNXq97tw4UJmzpyZpY6q5lRzEJk10U1zi3VLHWXxqZmZZecjIrM2rTp+8F7wZVi07Lm9J1YdX11XrEXLqj8+cRCZtWHV8cv55i6P5y6jEv9+xVqVh5FPzcwsOweRmWXnIDKz7BxEZpadg8jMsvNdMytd/cSGnuTQBuMgstLlnNzRxgafmplZdj4iMmvi6aeffvb7RcvEv1+xVsZqqlP/VHf9n0GZHERmTTx3hhs9r9tFP6hqlh+fmplZdj4iMmtCUt0RQVTa8TSndOSXjv6kao4CHURmTUycOPHZcYBWHR991uk1BdDEiRMr2adPzcwsOweRmWXnIDKz7PouiCTtL+lySY9LWiDpFklHSpowzO1sJeloSZdIekjSM8U2/yDpkOFuz6yf9dXFakknAUcAS4FLgfnAbsAJwF6S9oyIIZ/gkrQSMLt4OR+4HngImArsCOwMHCTpjRHxz47/IGY9pm+OiCTtQwqh+cAOEfHGiNgPeAlwKyk8jh3GJm8E3gasExG7RcQ7IuK1wFbAg8D2wJc7+TOY9aq+CSLgmGJ5fETUjmaIiEeAjxQvD5U0ZagNRcTSiNg2Is6NiMUN624Fjixevt2naGZD64sgkrQhsF3x8uzG9RFxFTAArAJ0Yia7m4rlRGCdDmzPrKf1yzWirYrlYxFxb5M2NwDTirY/GuX+XlIslwCPjXJbAEiaOkST9TuxH7Mc+iWINi2W97doM9DQdkSUnomvnZr9qvHUbRQGhm5iZVm0bFylve89r1lvmlwsF7RoM79YrjHKfX2GdOdsPnDUKLdlXSTHL+iKffd2z/9+CaJKSDoI+DSwHHhfRNzVwc1PG2L9+qTHCMzGnH4JoqeK5Wot2qxeLJ8cyQ4k7Q+cUbz8YEScO5LtNBMRc4fYfyd3Z+QdX7vW2bYmVy1V7bdfgmhOsWx1VFFbN6dFm0FJ2pd0N24c8KGIOGOIj9gYkHOs7ZkzZz4bRpMmTeKiiy7KVksV+uL2PStup68tqdnF6G2L5ewm6wdVPCj5Y2A88OGI+NbISjTrX30RRMVpTe36yTsb10vamXREtBho+78eSXsB55COLD8cEd8YfbVm/acvgqhwXLE8StLWtTclrQ2cXrw8NSKeqFs3S9Ltki5p3JikmcB5pBD6D4eQ2cj1yzUiIuICSScDhwPXFuGyANgdWBO4GvhUw8emANOBVevflPRC4HxgZWAuMEPSjCa7/njRjcTMmuibIAKIiCMkXQ0cAswAJgB3A8cDX4mIJW1uahKpOwikHvcHt2j7WcBBZNZCXwURQEScQ7qu007bM4EzB3l/DrXRxc1s1PrpGpGZdSkHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXnIDKz7BxEZpadg8jMsnMQmVl2DiIzy85BZGbZOYjMLDsHkZll5yAys+xWyl2AmT3fpEmTBv2+VykictdgHSBpKjAAMDAwwNSpUzNXZH1AndqQT83MLDsHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXZ9F0SS9pd0uaTHJS2QdIukIyVNGOH2tpF0rqSHJC2SdK+kUyS9sNO1m/WqvgoiSScB5wA7AdcBvwE2Ak4ALpU0cZjbeytwLfBW4D7g58By4FDg/yS9uHPVm/WuvgkiSfsARwDzgR0i4o0RsR/wEuBWYGfg2GFs70XA90ijXH4oIraPiAOAlwI/ANYDzpbUscGjzHpV3wQRcEyxPD4iZtfejIhHgI8ULw+VNKXN7X0UmARcHBHfrNveMuDDwBPAdsCeoy3crNf1RRBJ2pAUCgBnN66PiKtIw6yuAsxsc7OzWmxvPvCL4uW+wyrWrA/1RRABWxXLxyLi3iZtbmho25SkyUDt+s8NTZq1vT2zftcvs3hsWizvb9FmoKFtK5vUfd9sm8PZ3pCKwfFb2bD2zYMPPtiJXZq1NG3atKnAvIhYOtpt9UsQTS6WC1q0mV8s1xjG9lptczjba8fA0E2S7bffvkO7NGtpAJgGzB3thvrl1MzMuli/HBE9VSxXa9Fm9WL55DC2V9vmE6PcXjumDbF+ZWAL4GHgH8CyDu23E9YHri++3w6Yl7GWsWKs/Jl1pK5+CaI5xbLVL3Nt3ZwWbWruq/t+I9JzSKPZ3pAiop3D33s6sa9Oa3iUal6bP0tf67c/s345NbupWK4tqdnF422L5ewm658VEU8Cf2v43Ii3Z9bv+iKIiv9Naoe572xcL2ln0hHMYuCiNjf7sxbbWx3Yq3h5/rCKNetDfRFEheOK5VGStq69KWlt4PTi5akR8UTdulmSbpd0ySDbOwlYCOwh6YN1nxlfbG9NUvj9rrM/hlnv6ZsgiogLgJNJF5GvlfRrSeeRTrFeCVwNfKrhY1OA6cDmg2zv78B7SBeFvynpWkk/Bu4EDgQeAt4ZEVHOT2TWO/omiAAi4gjgAOAaYAapO8dc4Chgt4h4epjbOxfYgXT6tRmp28d44DTgVRHxtxYfN7OC/B+2la14Krz2QOa0Xr8D1An99mfWV0dEZtadfERkZtn5iMjMsnMQmVl2DiIzy85BZGbZOYjMLDsHkZll5yAys+wcRGaWnYPIzLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIrMuI2lzSSdKukrSHZJOrFu3g6R/lzQlZ42dtlLuAmzsk7RsFB+PiPC/w3DffDsAAA1PSURBVIKkg4GvA6sUbwWwTl2TScDXgCXAmZUWVyIfEVknaBRf/jdYkPQa4NukkDmSNGeeGppdATzBiinNe4L/EdioRcS4xi/gy8Ai4KvA1sBaxddWpOm6nwa+XLS15EhS8LwpIr4UEdc3NoiI5cDNwMuqLq5MPiS2jpP0XuCjwBsi4rKG1bcA/yXpF8DFkm6LiO9UXmR32gm4LiKuGqLdPFK49wz/b2Rl+Ahw9SAh9KyIuBy4CvhwVUWNAWsC97fRbiKwcsm1VMpBZGXYghXTJbfyADC95FrGkkeBjdto92LSUVHPcBBZGZYCr2yj3SuKtpZcC2wr6eXNGkjaCXg56WiyZziIrAzXAq+QdHizBpIOI4XVNZVV1f1OA8YDP5X06saVkrYEziDd0j+94tpKpYjIXYP1GEk7AleS/qO7BjgbuLdYvQnwLmBHYBnw+oj4Y4Yyu5KkLwH/RQqbu4HNSadhj5LulI0DToyIo7IVWQIHkZVC0v6kZ2Imk36pnrMamA98MCJ+UnVt3U7SB4HPAhs0rHoE+FxEnFZ5USVzEFlpJK0HfADYBZhavP0A6aG870TEg7lq63aSRHrmajPSUdAAcH1E9OQ1NQeRmWXni9Vmlp2frLbSSFoDeDcwA1gXuCQiTizWvZR04frKiFiUrciMJL1uNJ+PiCs7VUtuDiIrhaQ9SXfL1iJdnA7S9aGa6cAFwDuAcyovsDtczvMv5Lcr6KHf3575Qax7FM+7/IzUDeFrpIvTjXfHfgssBPamf4PoSkYeRD3FQWRlOAZYFdg/Is4HkPScIIqIJZJuBl6Vob6uEBGvz11Dt/DFaivDrsAttRBqYS7Pf1amb0j6tKS35K6jGziIrAzrAne20W4lYLWSa+lmnwX2qb2QtExSXw6J4iCyMjwBbNhGu82Ah0uupZstJ/Utq6mNWtl3HERWhtnANpI2atZA0itI14f+VFlV3edR0pApfc8Xq60M3wb2BH4kab+IeM7YOZLWKdqoWPara4C9JF0J/K14b2dJZ7Tx2YiI95dXWrXcxcNKUdwl2580NvXVwB7AXcAdwOuB1YEfRsSBuWrMrXjM4SLaGwytUUTE+KGbjQ0OIiuFpPHA50ljV09sWL0EOAU4KiJGMxXRmCdpErAdKYzOJA141tYF64j4XnmVVctBZKWStBbpdn59L/JLIqKfL1IPStJy4MyIeF/uWqrmIDLLpLgWdFVEnFG8Phi4qx8HivNdM+s4SYcVR0LW2nuAnetenwH0zAXo4XAQWRm+Cvxd0k8k/WsxyJc93zJgQt1rP0dk1kHnk36h9gcuBAYkfUHSS/KW1XUeBl7toPY1IiuJpBeQxiJ6L+nBxdo/tKuB7wLnRMSCTOV1BUlnkSYSuI80ucDrSQPl397GxyMidi+vumo5iKx0kl4FvI809tA6pFBaCJxLukvUMwN8DYekF5GGS9luBB/3c0RmIyFpJeAtpKOkN5Ke7F8eEX39hL+kjUnPEV0O/AY4oZ3PRcQVJZZVKQeRVU7SVOATwCH02P/so9HPzxH19f9EVh1JqwCzSEdDu7HiRslfshXVfTYlzffWd3xEZKWStAPpeZkDgCmku2lPAD8GvhsR1+WrzrqFj4is4yRtABwIHEwa5qI2eP5lpDtmP+3XmTtakXTQcNpHxPfLqqVqPiKyjpP0DOnUS6Rb098jHf3cl7WwLldcI2rnF1L02LU1HxFZGZaSHmo8IyIuyV3MGPJ9Bg+icaS7aluThta9gHR62zN8RGQdJ2mNiHgydx29RtJ6wFnAC4EZEbEwc0kd4y4eVobHJc3OXUSviYiHSE9ibwx8KnM5HeUgsjIsBP6au4heFBH/AK4n9ePrGQ4iK8NdpNMHK8cC2pslZcxwEFkZfgC8VtLmuQvpNZKmADsC/8xdSyf5rpmV4STgdcClko4GzvdzQ0NrNf0SMBnYktQ1Zl3SHbae4btm1nGS7iE967IxK25HP0ya0aNRRISPnGj7OaLas1kzIuLB8quqhoPIOq74hWpXTz2YNxqS5tA8iJYADwAXA6dFRE89R+RTMyvDprkLGIsiYpPcNeTiIyIzy853zcwsO5+aWWkkrUEat3oG6U7PJRFxYrHupcAmwJW+o/ZcxdhN25KeFVq1Wbte6n3vILJSSNoTOBtYixXDgDxQ12Q6qfPmO4BzKi+wS0k6HPgsaeymoTiIzJqRtCVpUPiVga8BVwA/aWj2W1JXkL1xEAEg6UDSM1iQZvK4DeiLzsMOIivDMaRTiv0j4nwASc8JoohYIulm0lRDlnyUdOT43l467WqHL1ZbGXYFbqmFUAtzgQ0qqGes2BK4tt9CCBxEVo51gTvbaLcSaaAvSxYBc3IXkYODyMrwBO31Dt+M1PXDkhuAvpyW20FkZZgNbNOqE6ekV5CuD/2psqq63xdJf27/lruQqvlitZXh28CewI8k7RcR8+pXSlqnaKNi2ZcGCeq7gf8GfibpZOBXwP3AoH33IuL+ciusjrt4WCmKu2T7k3rcXw3sQRow7Q7g9cDqwA8j4sBcNebWore9mrxfL3ppqm4HkZVC0njg86Rb0hMbVi8BTgGOiohlVdfWLYbobT+kiOiZzsUOIiuVpLVIt/M3I12THCB19fBFanuWg8jMsvNdM6ucpCmSlLuObidpJUkfk/QHSbdJ+r2k9+WuqwwOIus4Sa+QdHjRw77+/V0l3Qs8Bjws6T1ZCuwSkvaV9LCkLwyybhxwIXAisBOpk/DuwLcknVlpoRVwEFkZDge+TN0Y1ZLWJvW235h0V2ht4NuStspSYXfYlfTncN4g6z4IvIH0Z/UL4FBSKD0NHFiMbtAzHERWhp2Av0TEQN17B5JmovgGsCZwEOnf32HVl9c1dgAejIibBln3IdIdtR9FxD4RcXpEHEX6cxPpz7Nn+GK1dZykR4BrImKvuvcuJD3kuH5EPFq8dyMwKSK2zFNpXpIGSJ2D39zw/jqkri8BbBMRNzesvwdYGhHPOfUdy3xEZGVYg9TfrN4OwM21ECrcRY/NWDpM6wCPD/L+dsXyH40hVPgr8KLSqsrAQWRleJK6gCkGSnsB8MdB2vbzIfky0kgFjbYulrObfO6f9Fj3LAeRleFmYIakFxev308KnCsa2m0K9MwkgSNwH7C1pJUb3t+d9OfVrEPwOsBDZRZWNQeRleEbwATgRkk3Af9JuuZxYa2BpMnAq4E/Z6mwO1xGumt2bO0NSbsCuxQvLxzsQ8BWwN/LLa1aDiLruIg4lzQA/EqkoT7uIw0bu7iu2dtIYdV4lNRPTiL1u/u4pAFJs0ljeQP8KSJuaPyApB1Jp3M9NXyK75pZaYpTjjUi4pFB1m1EmuHj7oiYX3lxXULSvsCZpNEIah4AdouIuwZpfybpFv7eEfHLKmqsgoPILDNJLwTeDKxHGn/ogohY0KTtR0hHkt9u1mYschBZqYqjoq2BqaQH8eYCN0bEkqyFWVdxEFkpigD6DPAR0nNF9Z4izXf22YbrRtanHETWcZJWBX5H6uoh4FFWzE6xCelOUZCeK3qDp5w23zWzMnwC2Jn05PReEbFuRGxXfK1Luh5yJzADODJjndYlfERkHSfpDtLEiVtExKDPu0h6EWn86r9HxPQq67Pu4yMiK8NGwGXNQgigWHdp0db6nIPIyvA4dWMRtbCIwTt9Wp9xEFkZLgZeJ2mVZg2KC9qvJR0VWZ9zEFkZPkl66O7s4mG95yjG2/lB0eb/VVybdSFfrLZRk/TpQd7ejNQVYSHpVv69xfubkgZImwh8H7gnIo4d5PPWRxxENmp1M5bWz8zR+LqZiIjxpRRmY0ZPDa5k2XwudwE2tvmIyMyy8xGRlUbSxqwYCvUfEXFfznqse/mIyDpK0nRSt429SH3K6j0K/BL4UkTcVnVt1r0cRNYxkg4FvkS6Ld/sQnUAS4EjI+KrVdVm3c1BZB1RDNh1CimAbgHOAq4nDfIu4IXA9qSJAf+FFEhHRMSpWQq2ruIgslGTNI3Um348KVy+NkT7Q0jjNS8FpkfE/eVXad3MT1ZbJxwKrAJ8YqgQAoiI00hDhawCHFJybTYG+IjIRq2YMuhFwAYRsbzNz4wnTYnzYES8usz6rPv5iMg6YWPSXPdthRBARCwDrik+a33OQWSdMJHUp2y4FgKrdrgWG4McRNYJ/wA2H8HnNgeeN+eZ9R8HkXXCDcA2krZo9wOSXgZsS7rFb33OQWSd8BPSv6WzJDVOHfQ8RZuzipc/LrMwGxscRDZqEfET0pHN1sCNkvaW9Lx/W5LGSZoFzAZeDdwQEedUW611I9++t46QtAFwFWngswD+CdxEerIa0nTKWwNTSE9azwF2iogHKy/Wuo6DyDpG0prAacABrDjarv0Dq/U9Ww6cAxwSER443wAHkZVA0qak3vfbsGIYkEeAG4FfRsQ9uWqz7uQgMrPsfLHazLJzEJlZdg4iM8vOQWRm2TmIzCw7B5GZZecgMrPsHERmlp2DyMyycxCZWXYOIjPLzkFkZtk5iMwsOweRmWXnIDKz7BxEZpadg8jMsnMQmVl2/x8uWZonr01ovQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 200x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot accuracy\n",
    "shuff_v=pd.DataFrame()\n",
    "shuff_v['accuracy'] = out['acc']\n",
    "shuff_v['shuff_accuracy'] = out['sacc']\n",
    "fig = plt.figure(figsize=(1, 2))\n",
    "fig.dpi = 200\n",
    "ax = sns.boxplot(data=shuff_v, showfliers=False, showcaps=False, order=['accuracy', 'shuff_accuracy'])\n",
    "plt.xticks(ticks=[0, 1], labels=['Observed', 'Shuffle'], rotation=90)\n",
    "ax.set(ylabel='Accuracy')\n",
    "plt.ylim([0,1])\n",
    "plt.axhline(y=shuff_v['shuff_accuracy'].quantile(.95), xmin=0, xmax=1, ls='--', c='k') #plot threshold for significance \n",
    "sns.despine()\n",
    "fig.savefig(data_loc +decoder+'accuracy.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6b0891b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95th precentile of suffle is:  0.3275510204081632\n",
      "mean of decoding distribution is:  0.37693877551020405\n",
      "mean shuffle accuracy is:  0.23755102040816325\n"
     ]
    }
   ],
   "source": [
    "# is it significant?\n",
    "loc = int(np.floor(len(shuff_v['shuff_accuracy'])*0.95))\n",
    "a=np.argsort(shuff_v['shuff_accuracy'])\n",
    "print('95th precentile of suffle is: ',shuff_v['shuff_accuracy'].quantile(.95)) # [a].to_numpy()[loc])\n",
    "print('mean of decoding distribution is: ',np.mean(np.array(shuff_v['accuracy'])))\n",
    "print ('mean shuffle accuracy is: ',np.mean(shuff_v['shuff_accuracy']))"
   ]
  }
 ],
 "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.7.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}