Skip to content

Instantly share code, notes, and snippets.

@smsharma
Created November 2, 2023 03:46
Show Gist options
  • Save smsharma/d3aea2558e66c5361a829d04307ba38d to your computer and use it in GitHub Desktop.
Save smsharma/d3aea2558e66c5361a829d04307ba38d to your computer and use it in GitHub Desktop.
GP fit with background test
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "75ccbc23-d4a1-45b1-8f06-411f57a45c46",
"metadata": {},
"outputs": [],
"source": [
"import jax\n",
"import jax.numpy as np\n",
"import numpy as onp\n",
"from scipy.stats import poisson\n",
"from scipy.stats import chi2\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1586034b-8edf-40c8-8656-019e6e84848e",
"metadata": {},
"outputs": [],
"source": [
"def bump_forward_model(x, amp_s, mu_s, std_s, amp_b, exp_b):\n",
" \"\"\" Forward model for a Gaussian bump (amp_s, mu_s, std_s) on top of a power-law background (amp_b, exp_b).\n",
" \"\"\"\n",
" y_b = amp_b * (x ** exp_b) # Power-law background\n",
" y_s = amp_s * np.exp(-((x - mu_s) ** 2) / (2 * std_s ** 2)) # Gaussian signal\n",
"\n",
" y = y_b + y_s # Total mean signal\n",
"\n",
" return y"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d753e559-96f8-4897-b145-49c8b1f1ef9e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x137928490>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGwCAYAAACtlb+kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3uklEQVR4nO3deVhU1f8H8PcMO8giKgKymruVormgoiakpZmKlKUppuXXNdRcctdyyxY1I7dKs69bIm5pqLkk7iuuhBsqomhJgqDIMuf3h1/uj5FtBmbmzgzv1/PM88yce7j3M3Nn+XDuWRRCCAEiIiIiM6aUOwAiIiIifWPCQ0RERGaPCQ8RERGZPSY8REREZPaY8BAREZHZY8JDREREZo8JDxEREZk9S7kDMBYqlQp37tyBo6MjFAqF3OEQERGRBoQQePToETw9PaFUFt+Ow4Tnf+7cuQNvb2+5wyAiIqIySEpKgpeXV7HbmfD8j6OjI4BnL5iTk5PM0RAREZEm0tPT4e3tLf2OF4cJz//kX8ZycnJiwkNERGRiSuuOwk7LREREZPaY8BAREZHZY8JDREREZo99eIiowsjLy0NOTo7cYRCRFqysrGBhYVHu/TDhISKzJ4RASkoKHj58KHcoRFQGLi4ucHd3L9c8eUx4iMjs5Sc7bm5usLe35+SiRCZCCIHHjx/j/v37AAAPD48y74sJDxGZtby8PCnZqVKlitzhEJGW7OzsAAD379+Hm5tbmS9vsdMyEZm1/D479vb2MkdCRGWV//ktTx88JjxEVCHwMhaR6dLF55cJDxEREZk9JjxERERk9pjwEBERaUGhUGDz5s1yh0FaYsJDRGSE+vfvD4VCgcGDBxfaNmzYMCgUCvTv39/wgZmoipak3LhxAwqFAnFxcXKHYjSY8BARGSlvb2+sW7cOT548kcqysrKwZs0a+Pj4yBgZkelhwqNHmZmZUCgUUCgUyMzMlDscInpOZmZmsbesrCyN6xZMSIqrWxZNmjSBt7c3oqOjpbLo6Gj4+PggICBAra5KpcKcOXPg7+8POzs7NGrUCFFRUdL2vLw8DBw4UNpet25dLFy4UG0f/fv3R/fu3fHVV1/Bw8MDVapUwbBhw0odCrxlyxY0adIEtra2qFmzJmbMmIHc3FwAwGeffQZPT088ePBAqt+lSxe8+uqrUKlUAJ61vixevBhvvPEG7OzsULNmTbXYASApKQnvvPMOXFxc4Orqim7duuHGjRtqdX766Sc0bNgQNjY28PDwwPDhwwEAfn5+AIAePXpAoVBIj0uLHQCuXLmCtm3bwtbWFg0aNMDu3btLfC2AZ+di3rx5qFWrFmxsbODj44NZs2ZJ28+fP48OHTrAzs4OVapUwaBBg5CRkSFtb9++PUaOHKm2z+7du6u16Pn5+WH27NkYMGAAHB0d4ePjg2XLlknb/f39AQABAQFQKBRo3749AGD//v1o3rw5HBwc4OLigtatW+PmzZulPiezIEgIIURaWpoAINLS0nS2z4yMDAFAABAZGRk62y8Rae7Jkyfi0qVL4smTJ4W25X8+i7p17txZra69vX2xddu1a6dWt2rVqoXqaCs8PFx069ZNfPPNNyI4OFgqDw4OFvPnzxfdunUT4eHhUvnMmTNFvXr1RExMjLh27ZpYsWKFsLGxEfv37xdCCJGdnS2mTp0qTpw4Ia5fvy7++9//Cnt7e7F+/Xq1Yzo5OYnBgweL+Ph4sW3bNmFvby+WLVtWbJwHDhwQTk5OYuXKleLatWti165dws/PT0yfPl0IIURubq4IDAwU3bt3F0II8d133wkXFxdx8+ZNaR8ARJUqVcTy5ctFQkKCmDx5srCwsBCXLl2SYq9fv74YMGCAOHfunLh06ZLo3bu3qFu3rnj69KkQQojvv/9e2NraigULFoiEhARx/PhxMX/+fCGEEPfv3xcAxIoVK8Tdu3fF/fv3NYo9Ly9PvPjiiyI4OFjExcWJP//8UwQEBAgAYtOmTcW+JuPGjROVK1cWK1euFFevXhWxsbFi+fLlQohnvwseHh4iNDRUnD9/XuzZs0f4+/urnct27dqJiIgItX0+f759fX2Fq6uriIyMFFeuXBFz5swRSqVS/PXXX0IIIY4fPy4AiD/++EPcvXtXPHjwQOTk5AhnZ2cxZswYcfXqVXHp0iWxcuVKtXNhrEr6HGv6+82E53+Y8BCZJ1NPeO7fvy9sbGzEjRs3xI0bN4Stra34+++/1X4As7KyhL29vTh8+LDaPgYOHCjee++9Yo8xbNgw0bNnT7Vj+vr6itzcXKns7bffFr169Sp2H8HBwWL27NlqZb/88ovw8PCQHl+7dk04OjqK8ePHCzs7O7F69Wq1+gDE4MGD1cpatGghhgwZIu2vbt26QqVSSdufPn0q7OzsxM6dO4UQQnh6eopJkyYVG2dRSUppse/cuVNYWlqK5ORkafvvv/9eYsKTnp4ubGxspATnecuWLROVK1dW+03Yvn27UCqVIiUlRQihecLz/vvvS49VKpVwc3MTixcvFkIIkZiYKACIM2fOSHUePHggAEhJsCnRRcLDpSWIqMIqeBnhec9PX5+/lk9RlEr13gHPX2opj2rVqqFLly5YuXIlhBDo0qULqlatqlbn6tWrePz4MV577TW18uzsbLVLX5GRkfjpp59w69YtPHnyBNnZ2WjcuLHa3zRs2FDtuXt4eOD8+fPFxnf27FkcOnRI7ZJNXl4esrKy8PjxY9jb26NmzZr46quv8J///Ae9evVC7969C+0nMDCw0OP8Drdnz57F1atX4ejoqFYnKysL165dw/3793Hnzh0EBwcXG2dZYo+Pj4e3tzc8PT2LjfN58fHxePr0abGxxMfHo1GjRnBwcJDKWrduDZVKhYSEBFSvXl3j+F9++WXpvkKhgLu7e4nvU1dXV/Tv3x+dOnXCa6+9hpCQELzzzjvlWp/KlDDhIaIKq+CPjlx1NTFgwACpP0pkZGSh7fmJ2/bt21GjRg21bTY2NgCAdevWYcyYMfj6668RGBgIR0dHfPnllzh27JhafSsrK7XHCoVC6mtTlIyMDMyYMQOhoaGFttna2kr3Dxw4AAsLC9y4cQO5ubmwtNT85ycjIwNNmzbF6tWrC22rVq1aoYRTm/1qErs28td9Kg+lUgkhhFpZUf2otD1XALBixQp8/PHHiImJwfr16zF58mTs3r0bLVu2LHfcxo6dlomIjNzrr7+O7Oxs5OTkoFOnToW2N2jQADY2Nrh16xZq1aqldvP29gYAHDp0CK1atcLQoUMREBCAWrVq4dq1a+WOrUmTJkhISCh03Fq1akmJyPr16xEdHY39+/fj1q1b+Pzzzwvt5+jRo4Ue169fXzrGlStX4ObmVugYzs7OcHR0hJ+fH/bs2VNsnFZWVsjLy9Mq9vr16yMpKQl3794tNs7n1a5dG3Z2dsXGUr9+fZw9e1atI/uhQ4egVCpRt25dAM+SuILHzMvLw4ULF0o87vOsra2lv31eQEAAJkyYgMOHD+PFF1/EmjVrtNq3qWILDxGRkbOwsEB8fLx0/3mOjo4YM2YMRo0aBZVKhTZt2iAtLQ2HDh2Ck5MTwsPDUbt2baxatQo7d+6Ev78/fvnlF5w4cUIazVNWU6dOxZtvvgkfHx+EhYVBqVTi7NmzuHDhAmbOnInbt29jyJAh+OKLL9CmTRusWLECb775Jt544w21VoUNGzbglVdeQZs2bbB69WocP34cP/74IwCgT58++PLLL9GtWzd89tln8PLyws2bNxEdHY1x48bBy8sL06dPx+DBg+Hm5oY33ngDjx49wqFDhzBixAgAkBKi1q1bw8bGBpUrVy419pCQENSpUwfh4eH48ssvkZ6ejkmTJpX4etja2mL8+PEYN24crK2t0bp1a/z999+4ePEiBg4ciD59+mDatGkIDw/H9OnT8ffff2PEiBHo27evdDmrQ4cOGD16NLZv344XXngB33zzDR4+fKjVeXFzc4OdnR1iYmLg5eUFW1tbpKamYtmyZXjrrbfg6emJhIQEXLlyBf369dNq3yZLP92LTA87LROZp5I6Oxqz/E7LxXm+E6tKpRILFiwQdevWFVZWVqJatWqiU6dO4s8//xRCPOvY3L9/f+Hs7CxcXFzEkCFDxKeffioaNWpU4jEjIiIKdcp+XkxMjGjVqpWws7MTTk5Oonnz5mLZsmVCpVKJ4OBg0alTJ7UOxyNGjBAvvPCCePTokRDiWYfiyMhI8dprrwkbGxvh5+enNnpMCCHu3r0r+vXrJ6pWrSpsbGxEzZo1xUcffaT2nb1kyRLp+Xt4eIgRI0ZI27Zu3Spq1aolLC0tha+vb6mx50tISBBt2rQR1tbWok6dOiImJqbUUVp5eXli5syZwtfXV1hZWQkfHx+1ztHnzp0Tr776qrC1tRWurq7io48+kl4LIZ6NShsyZIhwdXUVbm5uYs6cOUV2Ws4fhZavUaNGYtq0adLj5cuXC29vb6FUKkW7du1ESkqK6N69u/Dw8BDW1tbC19dXTJ06VeTl5RX7XIyFLjotK4R47kJhBZWeng5nZ2ekpaXByclJJ/vMzMxEpUqVADy7Vqzr6/pEVLqsrCwkJibC39+/zP0ySL8UCgU2bdqE7t27yx0KGamSPsea/n6zDw8RERGZPSY8REREZPbYaZmIiGTFnhVkCGzhISIiIrPHhIeIiIjMHhMeIiIiMntMeIiIiMjsMeEhIiIis8eEh4iISAsKhQKbN2+WOwxZtG/fHiNHjiyxjp+fHxYsWGCQeLTBhIeIyAj1798fCoUCgwcPLrRt2LBhUCgU6N+/v+EDM1EVLUm5ceMGFAoF4uLi5A7FaDDhISIyUt7e3li3bh2ePHkilWVlZWHNmjXw8fGRMTIi08OEh4jISDVp0gTe3t6Ijo6WyqKjo+Hj44OAgAC1uiqVCnPmzIG/vz/s7OzQqFEjREVFSdvz8vIwcOBAaXvdunWxcOFCtX30798f3bt3x1dffQUPDw9UqVIFw4YNQ05OTolxbtmyBU2aNIGtrS1q1qyJGTNmIDc3FwDw2WefwdPTEw8ePJDqd+nSBa+++ipUKhWAZ60vixcvxhtvvAE7OzvUrFlTLXYASEpKwjvvvAMXFxe4urqiW7duuHHjhlqdn376CQ0bNoSNjQ08PDwwfPhwAM8usQBAjx49oFAopMelxQ4AV65cQdu2bWFra4sGDRpg9+7dJb4WwLNzMW/ePNSqVQs2Njbw8fHBrFmzpO3nz59Hhw4dYGdnhypVqmDQoEHIyMiQthd12ah79+5qLXp+fn6YPXs2BgwYAEdHR/j4+GDZsmXSdn9/fwBAQEAAFAoF2rdvDwDYv38/mjdvDgcHB7i4uKB169a4efNmqc+poNzcXAwfPhzOzs6oWrUqpkyZUuLkkT/88ANcXFywZ88eAMCjR4/Qp08fODg4wMPDA/Pnz9foUlm56WNVU1PE1dKJzFNJqyxnZGSIjIwMtZW8nz59KjIyMkRWVlaRdQuuLJ2dnS0yMjIK7buoutrKX7n8m2++EcHBwVJ5cHCwmD9/fqHVs2fOnCnq1asnYmJixLVr18SKFSuEjY2N2L9/vxTr1KlTxYkTJ8T169fFf//7X2Fvb6+2Knl4eLhwcnISgwcPFvHx8WLbtm3C3t5ebfXw5x04cEA4OTmJlStXimvXroldu3YJPz8/MX36dCGEELm5uSIwMFB0795dCCHEd999J1xcXMTNmzelfQAQVapUEcuXLxcJCQli8uTJwsLCQly6dEmKvX79+mLAgAHi3Llz4tKlS6J3796ibt264unTp0IIIb7//ntha2srFixYIBISEsTx48el1cTv378vAIgVK1aIu3fvivv372sUe15ennjxxRdFcHCwiIuLE3/++acICAgodbX0cePGicqVK4uVK1eKq1evitjYWLF8+XIhxLP3hoeHhwgNDRXnz58Xe/bsEf7+/mrnsl27diIiIkJtn0Wtlu7q6ioiIyPFlStXxJw5c4RSqRR//fWXEEKI48ePCwDijz/+EHfv3hUPHjwQOTk5wtnZWYwZM0ZcvXpVXLp0SaxcuVLtXJSmXbt2olKlSiIiIkL89ddf0vuo4Huk4EruX3zxhahSpYo4duyYtP3DDz8Uvr6+4o8//hDnz58XPXr0EI6OjoWec0G6WC2dCc//MOEhMk8lfVHmfz7zfwCFeJY4ABAffvihWl17e3sBQCQmJkpl8+fPFwBE79691epWrVpVABAXLlwoc9z5Cc/9+/eFjY2NuHHjhrhx44awtbUVf//9t9oPYFZWlrC3txeHDx9W28fAgQPFe++9V+wxhg0bJnr27Kl2TF9fX5GbmyuVvf3226JXr17F7iM4OFjMnj1breyXX34RHh4e0uNr164JR0dHMX78eGFnZydWr16tVh+AGDx4sFpZixYtxJAhQ6T91a1bt1BiamdnJ3bu3CmEEMLT01NMmjSp2DiLSlJKi33nzp3C0tJSJCcnS9t///33EhOe9PR0YWNjIyU4z1u2bJmoXLmy2m/C9u3bhVKpFCkpKUIIzROe999/X3qsUqmEm5ubWLx4sRBCiMTERAFAnDlzRqrz4MEDAUBKgsuiXbt2on79+mrnYvz48aJ+/fpqsc2fP1+MGzdOeHh4qH0O0tPThZWVldiwYYNU9vDhQ2Fvb6/3hEf2S1oHDhxA165d4enpWWynsvj4eLz11ltwdnaGg4MDmjVrhlu3bknbs7KyMGzYMFSpUgWVKlVCz549ce/ePQM+CyIi/ahWrRq6dOmClStXYsWKFejSpQuqVq2qVufq1at4/PgxXnvtNVSqVEm6rVq1CteuXZPqRUZGomnTpqhWrRoqVaqEZcuWqX2XAkDDhg1hYWEhPfbw8MD9+/eLje/s2bP47LPP1I770Ucf4e7du3j8+DEAoGbNmvjqq6/wxRdf4K233kLv3r0L7ScwMLDQ4/j4eOkYV69ehaOjo3QMV1dXZGVl4dq1a7h//z7u3LmD4OBgDV9VzWKPj4+Ht7c3PD09i43zefHx8Xj69GmxscTHx6NRo0ZwcHCQylq3bg2VSoWEhASt4n/55Zel+wqFAu7u7iWeK1dXV/Tv3x+dOnVC165dsXDhQty9e1erYwJAy5YtoVAopMeBgYG4cuUK8vLypLKvv/4ay5cvx8GDB9GwYUOp/Pr168jJyUHz5s2lMmdnZ9StW1frOLQl++KhmZmZaNSoEQYMGIDQ0NBC269du4Y2bdpg4MCBmDFjBpycnHDx4kXY2tpKdUaNGoXt27djw4YNcHZ2xvDhwxEaGopDhw4Z8qkQkYnJ7zdhb28vlY0dOxYjR46EpaX612P+D4mdnZ1UNmzYMHz00UdqCQIAqW9JwbrlMWDAAKk/SmRkZKHt+c9j+/btqFGjhto2GxsbAMC6deswZswYfP311wgMDISjoyO+/PJLHDt2TK2+lZWV2mOFQiH1tSlKRkYGZsyYUeT3d8Hv6QMHDsDCwgI3btxAbm5uode3JBkZGWjatClWr15daFu1atWgVJbtf3dNY9eGLs65Uqks1CemqH5U2p4rAFixYgU+/vhjxMTEYP369Zg8eTJ2796Nli1bljvugoKCgrB9+3b8+uuv+PTTT3W677KSPeF544038MYbbxS7fdKkSejcuTPmzZsnlb3wwgvS/bS0NPz4449Ys2YNOnToAODZCa1fvz6OHj2q85NIROaj4H/Z+aytrWFtba1RXSsrq0I/OsXVLY/XX38d2dnZUCgU6NSpU6HtDRo0gI2NDW7duoV27doVuY9Dhw6hVatWGDp0qFRWsPWnrJo0aYKEhATUqlWr2Drr169HdHQ09u/fj3feeQeff/45ZsyYoVbn6NGj6Nevn9rj/I7ZTZo0wfr16+Hm5gYnJ6cij+Hn54c9e/bg1VdfLXK7lZWVWguEJrHXr18fSUlJuHv3Ljw8PKS4SlK7dm3Y2dlhz549+PDDD4vc58qVK5GZmSm9Tw4dOgSlUim1clSrVk2t5SUvLw8XLlwo9rkVJf89/PxzBp51ZA4ICMCECRMQGBiINWvWaPVb+XySfPToUdSuXVst8W/evDmGDx+O119/HZaWlhgzZgyAZ619VlZWOHHihDTSMC0tDZcvX0bbtm01jqEsZL+kVRKVSoXt27ejTp066NSpE9zc3NCiRQu1y16nTp1CTk4OQkJCpLJ69erBx8cHR44cKXbfT58+RXp6utqNiMgYWVhYID4+HpcuXSrUmgQAjo6OGDNmDEaNGoWff/4Z165dw+nTp7Fo0SL8/PPPAJ79EJ88eRI7d+7E5cuXMWXKFJw4caLcsU2dOhWrVq3CjBkzcPHiRcTHx2PdunWYPHkyAOD27dsYMmQIvvjiC7Rp0wYrVqzA7NmzCyUOGzZswE8//YTLly9j2rRpOH78uNSq1adPH1StWhXdunVDbGwsEhMTsX//fnz88ce4ffs2AGD69On4+uuv8e233+LKlSvS88+XnxClpKTg33//1Sj2kJAQ1KlTB+Hh4Th79ixiY2MxadKkEl8PW1tbjB8/HuPGjZMuKR49ehQ//vij9FxsbW0RHh6OCxcuYN++fRgxYgT69u2L6tWrAwA6dOiA7du3Y/v27fjrr78wZMgQPHz4UKvz4ubmBjs7O8TExODevXtIS0tDYmIiJkyYgCNHjuDmzZvYtWsXrly5gvr162u171u3bmH06NFISEjA2rVrsWjRIkRERBSq16pVK+zYsQMzZsyQJiJ0dHREeHg4xo4di3379uHixYsYOHAglEql2mUyvdCgj5LB4LmOYHfv3hUAhL29vfjmm2/EmTNnxJw5c4RCoZA6Xa1evVpYW1sX2lezZs3EuHHjij3WtGnTpA6LBW/stExkXkrq7GjM8jstF+f5TqwqlUosWLBA1K1bV1hZWYlq1aqJTp06iT///FMI8axjc//+/YWzs7NwcXERQ4YMEZ9++qlo1KhRiceMiIgQ7dq1KzHWmJgY0apVK2FnZyecnJxE8+bNxbJly4RKpRLBwcGiU6dOap1cR4wYIV544QXx6NEjIcSz7/7IyEjx2muvCRsbG+Hn56c2ekyIZ78H/fr1E1WrVhU2NjaiZs2a4qOPPlL7zl6yZIn0/D08PMSIESOkbVu3bhW1atUSlpaWwtfXt9TY8yUkJIg2bdoIa2trUadOHRETE1PqKK28vDwxc+ZM4evrK6ysrISPj49a5+hz586JV199Vdja2gpXV1fx0UcfSa+FEM9GpQ0ZMkS4uroKNzc3MWfOnCI7LeePhMrXqFEjMW3aNOnx8uXLhbe3t1AqlaJdu3YiJSVFdO/eXXh4eAhra2vh6+srpk6dKo0mzO/ovG/fvmKfW7t27cTQoUPF4MGDhZOTk6hcubKYOHGi2vl9PrY///xTODg4iG+//VYI8azjcu/evYW9vb1wd3cX33zzjWjevLn49NNPiz2u2Y3Sev5NlJycLAAUGmXQtWtX8e677wohyp7wZGVlibS0NOmWlJTEhIfIDJlqwlORlJZAkGHs3btXuLi4iNTUVIMeNyMjQzg7O4sffvih2DpmMUqrJFWrVoWlpSUaNGigVl6/fn1pZIG7uzuys7MLNffdu3cP7u7uxe7bxsYGTk5OajddK3jt9MCBA0VeSyUiIjIGO3bswMSJE1G5cmW9HufMmTNYu3atdOm1T58+AIBu3brp9bhGnfBYW1ujWbNmhYbqXb58Gb6+vgCApk2bwsrKSprBEQASEhJw69atUocP6lN0dLRaota5c2f4+fmpzZhKRERkLL788kuMHTvWIMf66quv0KhRI4SEhCAzMxOxsbGFplvQNdlHaWVkZODq1avS48TERMTFxcHV1RU+Pj4YO3YsevXqhbZt2+LVV19FTEwMtm3bhv379wN4Nn5/4MCBGD16NFxdXeHk5IQRI0YgMDBQthFa0dHRCAsLKzSsMDk5GWFhYYiKiipyGCQRUUX0/HclmbeAgACcOnXK4MdVCJnfafv37y9yqF14eDhWrlwJ4Nn6KHPmzMHt27dRt25dzJgxQ63pKysrC5988gnWrl2Lp0+folOnTvj+++9LvKT1vPT0dDg7OyMtLa1cl7fy8vLg5+cnjRx4nkKhgJeXFxITE4scbUFEupWVlYXExET4+/uXeW4VIpJXSZ9jTX+/ZU94jIWuEp7iErjn7du3T1rMjYj0hwkPkenTRcJj1H14TJGm03SXZTpvIiIiKhsmPDqWPxunruoRERFR+THh0bGgoCB4eXkVO2OkQqGAt7c3goKCDBwZERFRxcWER8csLCywcOFCACgy6RFCYMGCBeywTGRiMjMzoVAooFAokJmZKXc4RKQlJjx6EBoaiqioKHh6eqqVK5VKjB49mkPSicho9O/fH927d5cet2/fHiNHjtTrMRUKhdqaiESGwIRHT0JDQ3Hp0iXp8W+//Ybs7Gx8/fXXMkZFRGVl6JnT+/fvL7UoWVlZwd/fH+PGjUNWVpZejxsdHY3PP/9cr8cwBTdu3IBCoUBcXJzcoZCOMOHRo4KXrdq3b1/sZSw2lRMZN7lmTn/99ddx9+5dXL9+HfPnz8fSpUsxbdo0vR7T1dUVjo6Oej0GkRyY8BhYYmIizp07J3cYRKSh/JnTk5OT1crzZ07XZ9JjY2MDd3d3eHt7o3v37ggJCcHu3bul7SqVCnPmzIG/vz/s7OzQqFEjREVFSdvz8vIwcOBAaXvdunWlPobFKXhJa//+/dI/YwVv/fv3l+pv2bIFTZo0ga2tLWrWrIkZM2YgNzdX2n7lyhW0bdsWtra2aNCggVr8xVGpVJg3bx5q1aoFGxsb+Pj4YNasWdL28+fPo0OHDrCzs0OVKlUwaNAgZGRkFPkc8nXv3l0tbj8/P8yePRsDBgyAo6MjfHx8sGzZMmm7v78/gGezAisUCmnetP3796N58+ZwcHCAi4sLWrdujZs3b5b6nEh+THgMqGPHjqhZsyY++OADuUMhIg3k5eUhIiKiyKUP8stGjhxpkIWBL1y4gMOHD8Pa2loqmzNnDlatWoUlS5bg4sWLGDVqFN5//338+eefAJ4lDl5eXtiwYQMuXbqEqVOnYuLEifj11181OmarVq1w9+5d6bZ3717Y2tqibdu2AIDY2Fj069cPERERuHTpEpYuXYqVK1dKyYlKpUJoaCisra1x7NgxLFmyBOPHjy/1uBMmTMDcuXMxZcoUXLp0CWvWrEH16tUBPGsR79SpEypXrowTJ05gw4YN+OOPPzB8+HCtXk8A+Prrr/HKK6/gzJkzGDp0KIYMGSKt3Xj8+HEAwB9//IG7d+8iOjoaubm56N69O9q1a4dz587hyJEjGDRoULGjcsnIlHU5d3Oj6fLy2sjIyBAABACRkZEhRo0aJQCISpUqlViPiHTnyZMn4tKlS+LJkyda/+2+ffukz2ZJt3379uk87vDwcGFhYSEcHByEjY2NACCUSqWIiooSQgiRlZUl7O3txeHDh9X+buDAgeK9994rdr/Dhg0TPXv2VDtOt27dpMft2rUTERERhf7un3/+ETVr1hRDhw6VyoKDg8Xs2bPV6v3yyy/Cw8NDCCHEzp07haWlpUhOTpa2//777wKA2LRpU5HxpaenCxsbG7F8+fIity9btkxUrlxZ7bty+/btQqlUipSUlGKfQ7du3UR4eLj02NfXV7z//vvSY5VKJdzc3MTixYuFEEIkJiYKAOLMmTNSnQcPHggAYv/+/UXGRvpT0udY099v2RcPrUg++eQTLFy4EBkZGbh586a04jsRGSe5Z05/9dVXsXjxYmRmZmL+/PmwtLREz549AQBXr17F48eP8dprr6n9TXZ2NgICAqTHkZGR+Omnn3Dr1i08efIE2dnZaNy4sVZx5OTkoGfPnvD19VW7JHb27FkcOnRI7XJTXl4esrKy8PjxY8THx8Pb21ttxGpgYGCJx4qPj8fTp08RHBxc7PZGjRrBwcFBKmvdujVUKhUSEhKkliBNvPzyy9J9hUIBd3d33L9/v9j6rq6u6N+/Pzp16oTXXnsNISEheOeddziRrIngJS0DqlGjBlq1agUAiImJkTkaIiqN3DOnOzg4oFatWmjUqBF++uknHDt2DD/++CMASH1Wtm/fjri4OOl26dIlqR/PunXrMGbMGAwcOBC7du1CXFwcPvjgA2RnZ2sVx5AhQ5CUlIQNGzbA0vL//0/OyMjAjBkz1I5//vx5XLlypczrltnZ2ZXp7wpSKpWFLkPm5OQUqmdlZaX2WKFQQKVSlbjvFStW4MiRI2jVqhXWr1+POnXq4OjRo+WOmfSPCY+Bvf766wCA33//XeZIiKg0xjRzulKpxMSJEzF58mQ8efIEDRo0gI2NDW7duoVatWqp3by9vQEAhw4dQqtWrTB06FAEBASgVq1auHbtmlbH/eabb/Drr79iy5YtqFKlitq2Jk2aICEhodDxa9WqBaVSifr16yMpKUmtBay05KB27dqws7PDnj17itxev359nD17Vm1E66FDh6BUKlG3bl0AQLVq1dSOmZeXhwsXLmj1vPP7ShXVPysgIAATJkzA4cOH8eKLL2LNmjVa7ZvkwYTHwPI7+23btk1tVAERGZ+SZk7Pf2zImdPffvttWFhYIDIyEo6OjhgzZgxGjRqFn3/+GdeuXcPp06exaNEi/PzzzwCeJQ8nT57Ezp07cfnyZUyZMgUnTpzQ+Hh//PEHxo0bhy+//BJVq1ZFSkoKUlJSkJaWBgCYOnUqVq1ahRkzZuDixYuIj4/HunXrMHnyZABASEgI6tSpg/DwcJw9exaxsbGYNGlSice0tbXF+PHjMW7cOKxatQrXrl3D0aNHpZatPn36wNbWFuHh4bhw4QL27duHESNGoG/fvtLlrA4dOmD79u3Yvn07/vrrLwwZMgQPHz7U6rV2c3ODnZ0dYmJicO/ePaSlpSExMRETJkzAkSNHcPPmTezatQtXrlxB/fr1tdo3yYMJj4EFBgZKzabLly+XOxwiKkVxM6d7eXkhKirKoDOnW1paYvjw4Zg3bx4yMzPx+eefY8qUKZgzZw7q16+P119/Hdu3b5eGVP/nP/9BaGgoevXqhRYtWuDBgwcYOnSoxsc7ePAg8vLyMHjwYHh4eEi3iIgIAECnTp3w22+/YdeuXWjWrBlatmyJ+fPnS/0TlUolNm3ahCdPnqB58+b48MMP1fr7FGfKlCn45JNPMHXqVNSvXx+9evWS+tbY29tj586dSE1NRbNmzRAWFobg4GB899130t8PGDAA4eHh6NevH9q1a4eaNWvi1Vdf1fh5A89e62+//RZLly6Fp6cnunXrBnt7e/z111/o2bMn6tSpg0GDBmHYsGH4z3/+o9W+SR4K8fyFzgoqPT0dzs7OSEtLg5OTk072mZmZiUqVKgF4dq07v5Ndv3794OLigk8++QS+vr7F1iOi8svKykJiYiL8/f3L3K8E+P/vCADYsWMHOnbsyDXxiAykpM+xpr/fHKUlg1WrVullv0yciPSnYHLTtm1bJjtEJoYJDxGRBhwcHIqcgJCITAP78Mjk4MGDCAsLw9y5c+UOhYiIyOwx4ZHJ5MmTsXHjRrWOdiXhAqNERERlx4RHJn379gUAPHjwQG2hPSLSD16OIjJduvj8MuGRSf/+/eHq6oqsrCxpkToi0r382XQfP34scyREVFb5n9/nZ8fWBjsty8TCwgIdO3bEunXrsHv3brnDITJbFhYWcHFxUZvHhatbE5kGIQQeP36M+/fvw8XFpVyjI5nwyOiNN97AunXrsHnzZrlDITJr7u7uAFDiwpBEZLxcXFykz3FZMeGRUevWrQEAV65ckTkSIvOmUCjg4eEBNze3IheRJCLjZWVlpZN5r5jwyOiFF16AnZ0dnjx5IncoRBWChYUFJwwkqqDYaVlmO3fuRHJystxhEBERmTW28MgsKCiI8+oQERHpGVt4iIiIyOwx4TECY8eOlTsEIiIis8aExwg8ffpUus/ZYImIiHSPCY8RmDhxonT/1KlTMkZCRERknpjw6JGDgwOEEBBCwMHBodh6BSdTWr9+vSFCIyIiqlCY8BiBvLw86f6aNWvULnERERFR+THhkVl0dDQaNGggPU5LS4Obmxuio6NljIqIiMi8MOGRUXR0NMLCwgpNPJieno6wsDAmPURERDrChEcmeXl5iIiIKHFU1siRI9UudxEREVHZMOGRSWxsLG7fvl3sdiEEkpKSEBsba8CoiIiIzBMTHpncvXtXp/WIiIioeEx4ZOLh4aHTekRERFQ82ROeAwcOoGvXrvD09IRCocDmzZuLrTt48GAoFAosWLBArTw1NRV9+vSBk5MTXFxcMHDgQGRkZOg38HIKCgqCl5cXFApFsXXs7OwQFBRkwKiIiIjMk+wJT2ZmJho1aoTIyMgS623atAlHjx6Fp6dnoW19+vTBxYsXsXv3bvz22284cOAABg0apK+QdcLCwgILFy4EgGKTnsqVK8PCwsKQYREREZkl2ROeN954AzNnzkSPHj2KrZOcnIwRI0Zg9erVsLKyUtsWHx+PmJgY/PDDD2jRogXatGmDRYsWYd26dbhz506x+3z69CnS09PVboYWGhqKqKioQkmcp6cnlEol7ty5g8uXLxs8LiIiInMje8JTGpVKhb59+2Ls2LFo2LBhoe1HjhyBi4sLXnnlFaksJCQESqUSx44dK3a/c+bMgbOzs3Tz9vbWS/ylCQ0NxaVLl6THO3bswK1bt9CpUycAwOrVq2WJi4iIyJwYfcLzxRdfwNLSEh9//HGR21NSUuDm5qZWZmlpCVdXV6SkpBS73wkTJiAtLU26JSUl6TRubRS8bNW2bVtYWFjg/fffBwAsXboUKpVKo/0UnLPnwIEDnMOHiIjof4w64Tl16hQWLlyIlStXlti5tyxsbGzg5OSkdjMmXbp0AQDcu3cPP/30U6n1n1+ionPnzvDz8+NszURUoWRmZkKhUEChUCAzM1PucMiIGHXCExsbi/v378PHxweWlpawtLTEzZs38cknn8DPzw/As5XG79+/r/Z3ubm5SE1NVVuF3NQ4OzvD398fSqUSV65cKbH1prglKpKTk7lEBREREYw84enbty/OnTuHuLg46ebp6YmxY8di586dAIDAwEA8fPgQp06dkv5u7969UKlUaNGihVyh68SePXuQmZmJFi1aFNt6U9ISFfllXKKCiIgqOku5A8jIyMDVq1elx4mJiYiLi4Orqyt8fHxQpUoVtfpWVlZwd3dH3bp1AQD169fH66+/jo8++ghLlixBTk4Ohg8fjnfffbfIIeymxN/fX2q9eT6hyW+9mT59usZLVLRv317PERMRERkn2Vt4Tp48iYCAAAQEBAAARo8ejYCAAEydOlXjfaxevRr16tVDcHAwOnfujDZt2mDZsmX6CtlgNGm9yZ/LpzRcooKIiCoy2Vt42rdvX+KK4c+7ceNGoTJXV1esWbNGh1EZB00WGE1NTdVoX1yigoiIKjLZW3ioeJq2yri6uhY7ik2hUMDb25tLVBARUYXGhMeIadoqExERAaDwEhX5jxcsWMAlKoiIqEJjwmPESltgNL/1ZtKkSUUuUeHl5YWoqCiEhoYaIlwiIiKjxYTHiBVcYPR5z7feFLVERWJiIpMdIiIiMOExevkLjDo6OqqVF9V6U9QSFURERMSExySEhobi4MGD0uNRo0ax9YaIiEgLTHhMxAsvvCDd/+yzz9h6Q0REpAUmPCZI1wupEhERmTsmPCYqKytL7hCIiIhMBhMeE9S6dWs4ODjg5s2bcodCRERkEpjwmKD4+HioVCp88skncodCRERkEpjwmKAPP/wQAHDixAnk5ubKHA0REZHxY8JjgmbMmIEqVarg1q1b2LJli9zhEBERGT0mPCbIzs4OgwcPBgDMnz9f5miIiIiMHxMeEzV06FBYWFjg0KFDWLVqldzhEBERGTUmPCbK09MTPj4+AJ5d4iIiIqLiMeExYbNnz0b79u2xefNmuUMhIjIZmZmZUCgUUCgUyMzMlDscMhBLuQMgwMHBAUIIrf/u3XffxbvvvquHiIiIiMwLW3iIiIjI7DHhMQPLli2Dh4cH/vOf/8gdChERkVFiwmMGtm/fjpSUFGzYsEHuUIiIZJWXlyfdP3DggNpjqtiY8JiBRYsWwcrKCtnZ2TrZHzv0EZEpio6ORoMGDaTHnTt3hp+fH6Kjo2WMiowFEx4z4OPjgyFDhsgdBhGRbKKjoxEWFobk5GS18uTkZISFhTHpISY85mLEiBFyh0BEJIu8vDxEREQUOdo1v2zkyJG8vFXBMeExE7Vq1UKHDh3kDoOIyOBiY2Nx+/btYrcLIZCUlITY2FgDRkXGhgmPGXnttdek+xcvXpQxEiIiw7l7965O65F5YsJjRoYNGybdr127toyREBEZjoeHh07rkXliwmNGHB0dIYSAEAKVK1eWOxwiIoMICgqCl5cXFApFkdsVCgW8vb0RFBRk4MjImDDhMWO5ublyh0BEpHcWFhZYuHAhABRKevIfL1iwABYWFgaPjYwHEx4zNWTIENjb2+O3336TOxQiIr0LDQ1FVFQUPD091cq9vLwQFRWF0NBQmSIjY6EQZVm10gylp6fD2dkZaWlpcHJykjuccvP09MTdu3fh5+eHxMRErf42MzMTlSpVAgBkZGTAwcFBHyESEelc/nc5AOzYsQMdO3Ys1LLD7zjzounvN1t4zNT3338PALhx4wZOnTolczRERIZRMLlp27YtL2ORhAmPmerevTvef/99AMDkyZOlci4bQUSmit9fVB5MeMzYjBkzYGlpiZiYGOzatUvucIiIiGTDhMeM1axZE7169QIAvP3221CpVDJHREREJA8mPGZu9OjRAJ516lqzZo3M0RAREcnDUu4ASL+aNGmCgQMHon379nj//fd53ZuIiCokJjwVwA8//CB3CCXiEFEiItI3XtKqYC5duiR3CERERAbHhKcCefPNN9G8eXO5wyAiIjI42ROeAwcOoGvXrvD09IRCocDmzZulbTk5ORg/fjxeeuklODg4wNPTE/369cOdO3fU9pGamoo+ffrAyckJLi4uGDhwIDIyMgz8TIwfJ+AiIqKKSvaEJzMzE40aNUJkZGShbY8fP8bp06cxZcoUnD59GtHR0UhISMBbb72lVq9Pnz64ePEidu/ejd9++w0HDhzAoEGDDPUUTMaKFSu4ijoRVXh5eXnS/QMHDqg9JjMmjAgAsWnTphLrHD9+XAAQN2/eFEIIcenSJQFAnDhxQqrz+++/C4VCIZKTk4vdT1ZWlkhLS5NuSUlJAoBIS0vTyXMxVt9++60AIACIO3fuFFknIyNDqpORkVHsvjStVxpd7YeIzJsm3xWl1dm4caOoUaOGVAeA8PLyEhs3btR3+KQnaWlpGv1+y97Co620tDQoFAq4uLgAAI4cOQIXFxe88sorUp2QkBAolUocO3as2P3MmTMHzs7O0s3b21vfoRuFvn37SveHDBkiYyRERIYVHR2NsLAwJCcnq5UnJycjLCwM0dHRMkVGhmBSCU9WVhbGjx+P9957T1oRNSUlBW5ubmr1LC0t4erqipSUlGL3NWHCBKSlpUm3pKQkvcZuLKysrKT7W7Zswfr162WMhojIMPLy8hAREQEhRKFt+WUjR47k5S0zZjIJT05ODt555x0IIbB48eJy78/GxgZOTk5qt4ro6NGjcodARKR3sbGxuH37drHbhRBISkpCbGysAaMiQzKJhCc/2bl58yZ2796tlpy4u7vj/v37avVzc3ORmpoKd3d3Q4dqUvbv34/58+cXKmeHPiIyN3fv3tVpPTI9Rp/w5Cc7V65cwR9//IEqVaqobQ8MDMTDhw9x6tQpqWzv3r1QqVRo0aKFocM1KQX7PeWLjo5GgwYNpMedO3eGn58fr20TkUnz8PDQaT0yPbInPBkZGYiLi0NcXBwAIDExEXFxcbh16xZycnIQFhaGkydPYvXq1cjLy0NKSgpSUlKQnZ0NAKhfvz5ef/11fPTRRzh+/DgOHTqE4cOH491334Wnp6eMz8x0/PDDDwgJCWGHPiIyW0FBQfDy8oJCoShyu0KhgLe3N4KCggwcGRmMAUaMlWjfvn1qwwPzb+Hh4SIxMbHIbQDEvn37pH08ePBAvPfee6JSpUrCyclJfPDBB+LRo0daxaHpsDZT9/yQzYMHD0qPXVxcin29FQqF8Pb2Frm5ucXuS1cxEREVpbzD0jdu3CgUCoVQKBSFvt8UCgWHppsoTX+/FUIU0WW9AkpPT4ezszPS0tLMugNzUQt1NmzYUOM1tvbt24f27dsXuy9dxURE9DxNvitKqxMdHY2PP/5YrSXb29sbCxYsQGhoqB6jJ33R9Pdb9ktaJL8DBw5onOSxQx8RmbLQ0FC1f/B27NiBxMREJjsVABMeQpUqVTBs2DCN6rJDHxGZuoLrCrZt25brDFYQTHgIAPDZZ5/B2tq62O3s0EdERKaMCQ8BeDY79ZgxY4rclj+qYcGCBfxPiIiMmoODA4QQEEKwPyCpYcJDklmzZqFZs2ZSh798Xl5eiIqK4jVuIiIyWUx4KpjSZlE+fvy42ugFdugjImPBWeCpPJjwVCCazqJc8LJVXl4eL2MRkew4CzyVFxOeCqKssyj37NkTqamphgiRiKhInAWedIEJTwWQl5eHiIgIFDXHZH7ZyJEji2wezs7OxqhRo/QeIxFRUcrz/UVUEBOeCiA2Nha3b98udrsQAklJSYiNjS1y+6pVqxATE6Ov8IiIilXe7y+ifEx4KgBNZ0cuqt5//vMfAED//v3x999/6zSufOyISETFKc/3F1FBTHgqAE1nRy6q3syZM9GgQQPcu3cPTZs2hUql0mls7IhIRCUpz/cXUUFMeCqAoKAgeHl5SRMIPq+kWZTt7Owwf/58AEBSUhKGDBmis7jYEZGISlOe7y+igpjwVAAWFhZYuHAhABT60tBkFuWOHTuie/fuqFu3LqZPn66TmNgRkYg0Ud7vL6J8THgqiNDQUERFRcHT01OtXNNZlKOjoxEfH6+zZmN2RCQiTZX3+4sIYMJToYSGhuLSpUvSY21mUVYoFNJ/UyqVCv/973/L1dmYHRGJSBvl+f4iApjwVDgFm33btm2rdTNwVlYWfHx80LdvX/j4+Ejl2nY2ZkdEItJWeb+/qGJjwkNasbW1lfrYpKWlqW3TprMxOyISEZEhMeEhreTl5RXZ0RjQrrMxOyISUUGZmZnSpfPMzEy5wyEzxISHtBIbG1tivxptOhuzIyIRERmKpdwBkGnRdWfj0NBQhISEwNnZGcCzjogdO3Zkyw4REekUW3hIK/robMyOiEREpG9MeEgr7GxMRKbOwcEBQggIIeDg4CB3OGQgTHhIK+xsTEREpogJD2mtpM7GY8eOxe7du2WKjIiIqGjstExlUlRnYyEEunTpAgBo0KABRowYIWeIREREErbwUJk939n49ddfh6+vLwBg+vTpuHnzplyhERGVC+cFMj9MeKiQsnboUyqVOHXqFBo3bozU1FT06NEDjx8/1mOkRETq7ty5g3bt2iEoKAhJSUlyh0NGhAkP6VSVKlWwefNmVK1aFWfOnEFoaChUKpXcYRFRBXD58mW0bt0aBw4cwMGDBxEYGIhz587JHRYZCa0SHmbLpAlfX19ERUVBqVRi586d6Nq1a7n3qUnzMpugiSquCxcuoE2bNrhx4wZq1aqFBg0aIDk5GW3atMHevXvlDo+MgFYJT7169TB16lRepqBStWvXDmFhYQCAuLg45ObmyhwREZkzb29veHp6okmTJjh06BAOHjyIdu3aAQCqVq0qc3RkDLRKeHbv3o2dO3eidu3aWLlypZ5CInOxdu1ajB07Fnfu3IGVlRVbXYhIb5ydnbFz507s27cPbm5uqFy5Mnbu3IkDBw7g5Zdfljs8MgJaJTytWrXCsWPHMGfOHEyZMgVNmzbVaJFIqpiUSiWmTZsmPVapVEhISJAxIiIyJ4sXL8b8+fOlx9WrV4eTk5P02MbGBo0bN5YeHzp0CMOHD2eLcwVVpk7L/fr1Q0JCArp06YI33ngDYWFhSExM1HVsZGZatmyJF198EadPn5Y7FCIycbNmzcLQoUMxevRoHD16tNT6jx49Qo8ePRAZGYlu3bohIyPDAFGSMSnXKK2OHTviww8/xKZNm9CgQQOMGzeObyIq1uXLl5Gbm4tu3brh33//lTscIjJhc+bMAfBszq8WLVqUWt/R0RHLli2Dra0tduzYgfbt2+PevXv6DpOMiFYJz5IlSzBw4EC8/PLLcHZ2RnBwMGJjYzF48GAsXLgQJ0+eRIMGDXDy5El9xUsmbMeOHXBxccHt27fRtWtXdn4nIq3k5eWpPf7+++8xbdq0Yhczfl737t2xb98+VK1aFadOnUJgYCBSU1P1ESoZIa0SnlmzZiEtLQ39+vXDvn378PDhQ5w6dQqRkZEYNGgQ9u7di8GDB6N///56CpdMWatWrXDgwAG4uLjg0KFDCAsLw5MnT+QOi4hMRFRUlHR/1apVGDJkiNb7aNmyJY4cOQJfX18kJiZi6dKlugyRjJjW8/BERUVhzJgxaNOmDezs7ArVGThwIOLj43UWIOlWWWdR1pWXXnoJ27Ztg42NDX7//Xe89NJLJjsxIef9ITIcIQSWLFkCAJgxYwb69u1b5n3VqlULM2bMAPDsysXzLUdknnQ+07KbmxsneaIStWnTBpMnTwYAXLt2DVu2bJE5IiIydiqVCm+//Tbq1KmDwYMHl3t/vXr1gq+vLzp37sy+pxWEzldLVygU0mRPRMWZPHkykpKSUKdOHfTo0UPucIjIyFlYWGDs2LEYM2aMxn12SmJra4tr166pLYJM5k32tbQOHDiArl27wtPTEwqFAps3b1bbLoTA1KlT4eHhATs7O4SEhODKlStqdVJTU9GnTx84OTnBxcUFAwcOZMZuApYuXYpPPvlE7jCIyIToItnJx2SnYpE94cnMzESjRo0QGRlZ5PZ58+bh22+/xZIlS3Ds2DE4ODigU6dOyMrKkur06dMHFy9exO7du/Hbb7/hwIEDGDRokKGeAulQRESE3CEQkQwK9qM5cOCA2uPly5dj3bp1yMnJ0flxhRA4evQo1q1bp/N9k5ERRgSA2LRpk/RYpVIJd3d38eWXX0plDx8+FDY2NmLt2rVCCCEuXbokAIgTJ05IdX7//XehUChEcnJyscfKysoSaWlp0i0pKUkAEGlpabp/YmYqIyNDABAAREZGRpnrCCFEy5YtpXrLli3T67F0xdDH04QxxkRUmo0bN4oaNWpI710AwsvLS2zcuFFkZmYKV1dXAUBs2bJF58feu3evACAqV64sMjMzpXJ+lkxHWlqaRr/fsrfwlCQxMREpKSkICQmRypydndGiRQscOXIEAHDkyBG4uLjglVdekeqEhIRAqVTi2LFjxe57zpw5cHZ2lm7e3t76eyJUqtjYWOkcDho0CMuXL5c5IiIyhOjoaISFhSE5OVmtPDk5GWFhYRg9ejRSU1Ph7++PLl266Pz4bdu2hZ+fH/7991+sX79e5/sn42HUCU9KSgqAZ+ujFFS9enVpW0pKCtzc3NS2W1pawtXVVapTlAkTJiAtLU26JSUl6Th60oalpSWOHz+OESNGAHiW9OTPpEpE5ikvLw8REREQQhTall/2448/AgBGjBihlz43FhYW0qiv77//Xuf7J+Nh1AmPPtnY2MDJyUntRvJSKBRYuHAhRo8eDQCYOHEiQkNDZY6KiPQlNjYWt2/fLna7EAK5ubmwtbXFgAED9BbHwIEDYWNjg5MnT+L48eN6Ow7Jy6gTHnd3dwAotN7JvXv3pG3u7u64f/++2vbc3FykpqZKdch0KBQKfPXVV2jTpg0AYNeuXZzUj8hM3b17V6N67dq1g7Ozs97iqFq1Kt555x0AbOUxZ0ad8Pj7+8Pd3R179uyRytLT03Hs2DEEBgYCAAIDA6UlLvLt3bsXKpVKowXlqOz0NWuzQqHAn3/+id69e+P8+fOyzAhNRPrn4eGhUb33339fz5EAQ4cOBQCsW7cO//zzj96PR4Yne8KTkZGBuLg4xMXFAXjWUTkuLg63bt2CQqHAyJEjMXPmTGzduhXnz59Hv3794Onpie7duwMA6tevj9dffx0fffQRjh8/jkOHDmH48OF499134enpKd8To3JRKpVYvXo1/P395Q6FiPQkKCgIXl5eJc6tY2tri/fee0/vsbRo0QIBAQFwd3fHtWvXShwmT6ZJ9oTn5MmTCAgIQEBAAABg9OjRCAgIwNSpUwEA48aNw4gRIzBo0CA0a9YMGRkZiImJga2trbSP1atXo169eggODkbnzp3Rpk0bLFu2TJbnQ/qR/98XAJNde4uI1FlYWGDhwoUACk8omL9O3U8//WSQCQIVCgW2bt2Ka9euITk5GQ0aNJC2de7cGX5+foiOjtZ7HKRH+h8hbxo0HcdP2tHVXBadO3eW9jNw4ECRk5Ojt2Npyhjn6dBVTMb43Mh8FTUPj7e3t9i4caMssSgUCrVYAAiFQiEUCoUsMVHJzGIeHqJ8a9eule7/+OOP6NGjR6HOzGyCJjJNoaGhuHTpkvQ4v4uCoUdpajJMfuTIkfxuMVFMeMjoRUdHqzUvA8Bvv/0GX19fXLx4scg6bIImMi0FL1stX74czZo1M3hiockw+aSkJMTGxhowKtIVJjykV+VtdSluFlYAePDgARo1aoSVK1eWOFMrkx4i09OvXz+DL+6p6TB5TeuRcWHCQ3pT3laXkpqX81laWmLKlClsgiYyIxYWFmoDFQxF02HymtYj48KEh/SitPVxNEl6SmteBoCnT5+yCZrIzISGhqJGjRoGP25pw+QVCgW8vb0RFBRk4MhIF5jwkM7pquOfLpuN2QRNZNwKzqgvR+sOoD5M/nn5SdCCBQsMfqmNdIMJD+mcrjr+6bLZmE3QRMbtv//9r3S/WbNmssURGhqKqKioQksTeXl5ISoqiuv7mTAmPKRzuur4p0nzspeXl2xN0Jp0yM7MzJQmUOOaYETFS01NlTsESWhoKBISEqTHffv2RWJiIpMdE8eEh3ROVx3/SpuFFQAWLlxYbB3gWWuSPpqgOQyeSLdmzZoldwhqCn5nLF68mJexzAATHtI5XXb8y29efn5dtILNy8XVyXfnzh3tn0QJdNEhm4iIDIsJD+mcJi0z2rS6PD8L644dOwo1Lz9fJzIyEo6OjlAqlTq9nMWZWIl0S6VS4erVq3KHUaqSpscg08CEh/RCk5YZbRRMjtq2bVtkslSwLDw8HHfu3EFMTAwaNWoklWdlZWl13OdxJlYi3Tpy5Ahq166Nzp07yx1Ksfr27YvWrVvLHQaVExMe0htNWmb0qVKlSnjttdekx9OnT4eLiwt2795d5n1yJlYi3VqzZg0AyDLvjqa2bt2KI0eOmERLFBWPCQ/plSYtM4aQnZ2NuXPn4unTp3jzzTexY8eOMu3H2Gdi5QKqZEpycnLw66+/AgDeffddmaMpXtu2bQEAmzZtkjkSKg8mPFQhWFtb4+jRo/Dw8EB2djbefPNNTJ8+Hbm5uVrtx5hnYuXIMTI1f/zxB/755x+4ubmhXbt2codTrK5duwJgwmPqmPBQhdG4cWPcuHEDgwcPhhACM2bMQK1atXDs2DGN96HrDtm6wpFjZIryL2f16tULlpaWMkdTvDfffBPAs/5GvFxtupjwUIVibW2NxYsXY/Xq1bC1tcXNmzcRGBiI69eva7wPXXfILi+OHCNT9PjxY6nFpHfv3jJHUzJPT0+0aNECALBlyxaZo6GyYsJDFVLv3r2xfft22Nvbo3nz5qhZs6ZWfy93h+yCOHKMTNG2bduQmZkJf39/KZkwZvmfbV7WMl3G24ZIpGcdOnTAvXv31C5NxcXF4fbt21ITdkmMpUO2OYwcy8zMRKVKlQAAGRkZcHBwkDki0rdu3bph48aNyM7OLrZPnDHp0aMHoqKi1EZ+kmlhwkMVWv6PLPBsJFeHDh3w77//YtiwYVi0aJFJfBEb+8gxoqLY2tqa1NpUtWvXxvHjx+UOg8qBl7SI/iclJUVKcCIjI6FUKk1i0U99jBzj8HYyNAcHBwghIIRgCx/pBRMeov/x8fHBvXv38PHHH0Op/P+PxsaNG2WMqnS6HjnG4e2kbwMGDMD06dNx//59uUPRWlpaGtauXYvHjx/LHQppiQkPUQGWlpZYuHAh/vjjD6ksPDwcL774IlJTU2WMrGS6GjnG4e2kb7dv38bKlSsxY8aMci/1IofmzZujd+/e5ZqxneTBhIcqnMzMTCgUihIvVzVv3lztcWJiImxtbQ0RXpmVd+SYvoa3a/J6k2HJeU7Wr18PIQSCgoLg4+Nj0GPrwhtvvAGAo7VMERMeolLMmjULa9euhb29PYBnC5Aaa1N8eUaOcXg7GcLq1asBGP/cO8Xp0aMHgGfD6rWdqZ3kxYSHqBQRERF46623pMc9evRAjRo1sHjxYhmj0j1zGN5Oxi0+Ph5nzpyBpaUlwsLC5A6nTNq0aYNq1aohNTUVBw4ckDsc0gITHiItZGRkYM+ePcjNzcXYsWPlDqdMirucweHtpA8F328///wzAKBTp06oWrWqzJGVjYWFhfQPEPu0mRYmPERaqFSpEq5evYqWLVuqlW/btk2miHTHmBdGJfOQvzK6KVzOKmmYfP5lrc2bN0OlUskRHpUBEx4iLfn4+ODIkSNq/92999578PX1xYULF9TqmtJ8Nsa6MCqZj1atWsHDw0PtErEpCg4ORqVKlZCcnIyLFy/KHQ5piAkPmQRNJiUz9MRlHTt2VHt869Yttbk5THE+G2NbGLUgjvYyfUuXLkVSUpLaDOemyNbWFr/++itu3bqFl156Se5wSENcWoJIB3744QecO3dOGs4eHR2Nnj17FqqXP5+N3MlDSUJDQxESEgJnZ2cAz4a3d+zYkS07pBPm8j7KH55OpoMtPEQ68O6770qXg/Ly8vDhhx8WWa8889kYkrEsjErmpag5nogMhQkP6VVFXB8nNjYW//77b7HbOZ8NVVSfffaZ3CHoVExMDN544w189dVXcodCGmDCQ6RjnM+G6P+lp6dL95/v92bqkpKSEBMTg/Xr18sdCmmACQ+RjnE+G+PGzs+a0dXrFBUVJd1/fjoHU/fWW29BoVDg5MmTSEpKkjscKgUTHiIdK20+GwCcz4YqjPzJBoHC0x2YuurVq6NNmzYAns3JQ8aNCQ+RjpU0n01+Wf58NlOnTsXEiRM5eRmZpXPnzuHUqVNyh6FX+ZMQbty4UeZIqDRMeIjKqKQO2cXNZ+Pt7S0NSf/nn38we/ZszJkzBxYWFrzEQiantMteP/74owxRGVb+9BOxsbG4d++ezNFQSZjwEOlJaGgoLl26JD3esWMHEhMTpfl3VCqVNG9PviFDhuDWrVsGjdMYmNKM1KQZIQT2798vdxh65+Pjg2bNmkGlUmHTpk1yh0MlMPqEJy8vD1OmTIG/vz/s7Ozwwgsv4PPPP1ebz0EIgalTp8LDwwN2dnYICQnBlStXZIya6JmS5rNxc3PD4cOHcebMGansl19+gb+/P1555RVcvnzZoLHKxRRnpKbSKRQKnD59ukL0benVqxdeffVV1KhRQ+5QqARGn/B88cUXWLx4Mb777jvEx8fjiy++wLx587Bo0SKpzrx58/Dtt99iyZIlOHbsGBwcHNCpUydkZWXJGDmRZmrXri3db9SoEVQqFU6dOoWzZ8/KGJVhREdHIywsDMnJyWrl+TNSM+kxbRYWFggJCZE7DL0bPXo09u7di65du8odCpXA6BOew4cPo1u3bujSpQv8/PwQFhaGjh074vjx4wCete4sWLAAkydPRrdu3fDyyy9j1apVuHPnToX4z4L0Q65LLLGxsZg5cyZeffVVvP3221L50qVL1eYzMQd5eXmIiIgocvZdU5mRmoqWkZGB3NxcucMwGHMbfWaujD7hadWqFfbs2SM17589exYHDx6U1jFJTExESkqK2n8Rzs7OaNGiBY4cOVLsfp8+fYr09HS1GxEg7yUWpVKJSZMmYe/evVLZyZMnMXjwYLi6upY6L4opzWwdGxuL27dvF7udM1Kbrnnz5sHb2xsrV66UOxSDSklJwdatW8v0t5wfSv+MPuH59NNP8e6776JevXqwsrJCQEAARo4ciT59+gB49gYDns2HUFD16tWlbUWZM2cOnJ2dpZu3t7f+ngQZlZJab4zxEsuhQ4egVCrV4ly/fr3JtHwU93pzRmr56LMFMy8vDytWrEBKSgrs7Ox0tl9jl5SUBE9PT/Ts2ROpqalyh0NFMPqE59dff8Xq1auxZs0anD59Gj///DO++uortcmsymLChAlIS0uTbpwls2IoqfXGWC+xRERE4O7duxg5cqRUNnDgQDg5OWHIkCHIzs42aDzaKOn15ozU8tB3C+auXbtw+/ZtuLq6onv37jrZpynw9vbGiy++iNzc3DK38pCeCSPn5eUlvvvuO7Wyzz//XNStW1cIIcS1a9cEAHHmzBm1Om3bthUff/yxxsdJS0sTAERaWlq5YybjtHHjRqFQKAQAtZtCoRAKhULMmDGj0Laibvv27dP4mBkZGdLfZWRklLnO8/Xs7Oyk+9evX9fqddBlTCUp7fX+9ddfhZeXV5F18ut5e3uL3NxcncWky/0Y6/FKUto52bhxo1S3rO+T0NBQAUBERERovB9z8dlnnwkAokuXLlr/bUV6nXRN099vo2/hefz4MZRK9TAtLCykmWn9/f3h7u6OPXv2SNvT09Nx7NgxBAYGGjRWMl6atN7kz45cGmO4xBIXF4fXXnsNXbt2hb+/v1Q+c+bMUkcnGqJDtiav9yeffIL58+cDKNzpM/9x/ozUVH6GaMG8f/++1LoxcOBA6bj5zH2OpbCwMADPWrnS0tJkjoYK0XvqVU7h4eGiRo0a4rfffhOJiYkiOjpaVK1aVYwbN06qM3fuXOHi4iK2bNkizp07J7p16yb8/f3FkydPND4OW3jM2759+zRqvdHkZgwtPEXVW758uQAgrKysRGRkpMjKyipUZ+PGjaJGjRpqz8fLy0vtP3ttYiqOpq/3vn37iozJ29tb5zHpej/GerziaHNONI37+TpffvmlACCaN28uhND8/WZOGjRoIACIX375Rau/M5b3iSnS9Pfb6BOe9PR0ERERIXx8fIStra2oWbOmmDRpknj69KlUR6VSiSlTpojq1asLGxsbERwcLBISErQ6DhMe87ZmzRqNvuxdXV01vsSiCUMmPGPHjlWL3d3dXcydO1f8/fffQgjdX84oiaav95o1a4QQ///5AyB27NhR5GusSR1NVNSER9tzou1799GjR6Ju3boCgFi2bJlW7zdzMnXqVAFAdOvWTau/M5b3iSkym4THUJjwmDdN/7udMWOG9IWsiy9pQyY8Qghx7949MWPGDLX/qhUKhWjWrJnw9PQs9nmXJZkria5bE3TZUlBREx5DtPCcPn1afPzxxyI1NVV4eXkZ7P1mTM6dOycAiEqVKml1lcFY3iemiAmPlpjwmLfc3FyNO8hqeolFE4ZOePI9ffpUrFixQri5uQng2WUubX7sykub17u056brloKKmvDo8pzkK66OtsmVOVGpVGLdunXiwYMHWv2dsbxPTJHZdFom0gULCwupU3JpHWRLW/TTFFhbW6N///5ITk7G5MmTNR4erKsO2dq83iUx1qkCTJGuzokmKvIcSwqFAr169YKrq6vcodBzmPBQhREaGoqoqCh4enqqlXt5eSEqKkotoSlp0U9TYmlpic8//xxDhw7VqL4u57zR5vUuDmdj1i1dnJPiDBo0CCdPngSg+fuIcyyRITHhoQrFHFpvyiIoKAheXl4lrvljZ2cHNzc3nR63vK+3PloKNB0mba5T/Wt6TrQdTr5mzRrs378fQOnvN4VCAW9vbwQFBZXxWRiP4t4nq1atQuvWrbFx40YZo6OCmPBQhWMurTfaKOlyRr4nT56gYcOGuHr1qs6PnU/b11vXLQVyrpNmTEo7J2V5nSwsLNCvXz/pvqEunxmrixcv4vDhw9iwYYPcodD/MOEhqiBKupzRs2dP2Nvbo3r16qhVq5a0bf369bIuXaHLlgJ9rJNmjq1AZX2d3nzzTbUWQn1ePjMF+ZMQ/vbbb3jy5InM0RDAhIdIrzRZvdyQM9EWdTnjxo0biIqKwqNHj3DmzBlp28WLF/Huu+/C3t4eEydOLHExXn1h52fD0vZ1Kngp8ZVXXin0+lXUS8jAs9fDx8cHmZmZ2Llzp9zhEJjwEMlKjkssxV3OUCqVapeGtm/fDoVCgby8PMyZMwfe3t545513sG7dOmlpF0Ng52fD0eZ1io6ORuPGjaVtU6ZMKfK9WxEvIQPPEvL8Vp6oqCiZoyGACQ+RbPRxiUWXxo0bh9TUVHz55ZcIDAxEbm4uNmzYgPfeew92dnbYtGmTwWIxxs7P5kjT579lyxaEhYUVuoxnLO9dY5Gf8GzduhVPnz6VORpiwkMkA1O5xOLi4oIxY8bg8OHDiIuLk77As7Oz0aJFC6nelStX9N7qY0ydn82Vps9/9erVRv/eNQYtWrRAjRo18OjRI+zevVvucCo8JjxEMjDFSyyNGjXChg0bkJycjG+//VbtElPTpk1hZ2eHESNG4N69ezJGWbSKNEy6PDR5napVq4a///672H0Y43tXLkqlEn369EG3bt1QuXJlucOp8JjwEMmgLJdYDNm5uSSenp4YMWKE9PjixYt49OgRsrOz8d1338HLyws9evTA2rVrkZWVJUuMz6tIw6TLM3JMk9epT58+Gu2rol8ezPfFF19g8+bNaN26dYn1jOXzbc6Y8BDJQNtLLMY8f0zDhg2RlJSEESNGoEWLFsjNzcXmzZvRu3dvODg44JNPPpE7RAAcJq2p0l6natWqabSfin55UBvG/Pk2K/pZysv0cPHQikOXC3qWlbaLmZrS4pnnz58XERERQqlUCgBi0qRJ0vEuXLggrl+/XqaYdBV3/mcdgNixY0exK3brKiZNj6cr+ow7JydH1KtXr8QFQcuyCKkp0/T5Xb16VWzfvr1Qua4/3xURV0vXEhOeisMYEh4h/v+L7vkvu4JfdPmJkaY/Lpow1A/Qo0ePxMSJE8U///wjHa99+/YCgGjYsKHYsmWLePr0qcYx6SpuXa5OX1qdjRs3iho1aqidMy8vL73+iOnqtSyqzvLlywUAUalSpVLfu9ocy5Rp8vwKrh6fkpIilevj810RcbV0IiOnySUWU+zcnK9SpUqYNWsWbG1tpbL4+HgAz/r9dOvWDZ6enhgyZAh++eUXucLUG2OfdkBbmZmZmDp1KgDg888/5+VBLTRt2lS6v379eum+KX++TRETHiIZlTa/jDHPH1OWzrFXr17FunXrMHz4cLi7u+PBgwdYsmQJhgwZAkdHR6hUqmJnpDYlpjLtgDYWLFiAu3fvws/PD0OGDKnQsyhrS6n8/5/apUuXSu8BY/58myMmPEQyK2l+GXObP0ahUKBXr15YtGgRkpKSEBMTgyZNmgAAfH19pZFAKpUKffr0wZEjR+QMt8zM7T/3v//+G1988QUAYPbs2bCxsQFQcWdRLo9Lly5J593cPt/GjgkPURE0WQPLEMx5/hhLS0t06tQJp06dQnJyMtauXStt27BhA9asWYNWrVohMDAQkZGRRjm/T3HM7T/39evX49GjR2jatCl69eoldzgm77vvvgNg3p9vY8SEh8iIVZT5Yzw9PfHiiy9Kjx89egRnZ2cAwNGjRzF8+HC1Vdxv3rxp8Bi1YW7/uQ8bNgzbtm3Dt99+q3Z5hspm06ZNSE5OrjCfb2PBdy5VOMbSeqMpU58/piyv94cffoiHDx/ixo0b+Oabb9C8eXO1/jD79u2T7htyIVNNGft/7tqeE4VCgTfffBOtWrUyQHTmrXXr1nBwcMD58+cBmP7n25Qw4SEyARW1g6ivry9GjRqFY8eO4eDBg1J5wdl+W7ZsCRcXFwwfPhx37tzR2bHLM/Mt/3OvOLR9nyxZsgTJycl4/fXXpbKK+vk2NCY8RCaioncQbdy4sXTfysoKwLPWnTNnziAtLQ2RkZGoUaMGWrdujVmzZmH//v1lPpYuZr419f/cC/5w79q1S8ZIjFdZ3if+/v5FtqpV9M+3ITDhIapgTO2SXkmUSiX279+Pbt26Sau3Hz58GJMnT8arr76KGjVqaL1PXc6fY6r/uatUKgwdOlR6XL9+/TLvy5zebwWV930ihMCJEyf0GSI9hwkPEZm01q1bY/PmzTh69CguX76sts3b21u6n52djWbNmmHBggXF7kvb+XM0uZxhav+5CyEwfPhwrF69Wior+DpS+edZysnJQZMmTdC8eXPExcXpNLbyLB5r7pjwEJFJ0KSloODlo3PnzuHHH3+UHi9btgwnT57E5MmTpbINGzbgn3/+kR5rM3+OOS74KITAmDFjsHjxYrlDMWrlnWfJysoKdevWBQBERkbqJUYqjAkPEZmlmjVromHDhtLjGjVqoE6dOmpLXXzwwQeoVq0aXFxcsHLlSo3nxdmyZYtZLRuRb9q0afjmm28AAN9//73M0RgvXcyzNHz4cADA6tWrkZqaqpO4qGRMeIioQujRowcSEhLUJjCsXbs2ACAtLQ1VqlTReF6c1atXm9WyEcCzyyxHjx4FACxatAj9+vWTOSLjpYt5llq3bo2XX34ZT548wYoVK3QVGpWACQ8RVSgF+9CcOXMGf/75J/r3748uXbpI8+cUR6FQoFq1avj777+LrWNqy0bks7KywtatWxEVFSW1PlDRdDHPkkKhkF7n77//3ijnkzI3THiIqEJr27YtVqxYAaVSqTZ/TlGKatUpjlzLRmg7L8zJkyel52Vra4uePXvqNT5zoKt5lnr37g0XFxdcv34du3fv1k+wJGHCQ0QVSmmdn0NDQ7Fx48ZC8+fke/TokUbH0WbZCE1G1mhSR9uO1P/973/RvHlzjBkzplAyZ67DyXVFF/MsOTg44IMPPgAAbN26tdT6uhyBVRFHczHhISJ6TmhoKOLj46XHO3bswPXr1zF27FiMGjWqxMsZAODu7o42bdoYIlSJtvPCREVFITw8HEIIPH361JChmg1N51kqqdUtIiICu3btwqJFi/QfcEUnSAghRFpamgAg0tLS5A6FqEgZGRkCgAAgMjIy5A7H6OIRQrcxlbSvjRs3CoVCIW0v6mZlZSXef/99jWIqb53c3Fzh5eVVbCwKhUJ4e3uLnJwcsXfvXtG9e3ehVCoFADFgwACRl5dXrteqIivt3G3cuFHUqFFD7Xx4eXmJjRs3arUfXdbRpp4p0PT3my08RERayr+c8fxMzs7OzmjYsCEUCgVycnKQk5Ojtr1ly5YYO3YsTp48WWhbeWg6L0zt2rXRoUMHbN68GSqVCgMHDsSyZcu4Arqe6HLWbio/vsuJqEzKs7imOSjqcsaDBw9w4cIFXLt2DZ9++inGjx+v9jcXLlzAV199hWbNmsHJyQlBQUH4z3/+g02bNpUrFk07SD948AD29vYYMmQILl68iB9++MHoZ342VdrOxsxLWvpnKXcARGR6oqOj8fHHH0uPO3fuDC8vLyxcuNDo14nSpeKWjfD398ecOXMAQK1DaJ06dSCEwD///IN///0XBw8exMGDB7Fs2TK1/R47dgy7d+9GzZo18cILL6BmzZqws7MrNg5NO0jPmjULffv2hYuLi6ZPkcpIm9mY27dvrzYsvagkicqPLTxEpBU205fd6dOncfnyZfzzzz/466+/8NVXX8HSsvD/nbNmzcKUKVPQp08ftGzZEm5ubmpJTcEWnc8++wwdO3Ys8bj588IMHTqUyY6BaDsbc9++faWyguuYke4w4SEyEcYwTLi8iybSM0qlEnXr1sUnn3yCzMxMHDt2TG17fsvACy+8IE2EmJGRIW0vuP6XhYVFif2BtJkXhnRH29mYvb298cUXXwAAJk6ciPv37+sttoqKCQ8Raay8iyZSYdbW1mprfgHPZt6dOXMmYmJikJSUhCdPnuDbb7+Vtru7u0v3O3bsiISEBPz777/YsGFDoY7U2swLQ7pTltmYR48ejcaNGyM1NRUjR440UKQVBxMeItKYLhZNpNK1bNkSkyZNQq1atQA8mwE5/z4AVKtWTbr//vvv46WXXsKFCxcQFham0bwwpH9lmY3Z0tISy5cvh1KpxNq1a/H7778bLuAKwCQSnuTkZLz//vuoUqUK7Ozs8NJLL+HkyZPSdiEEpk6dCg8PD9jZ2SEkJARXrlyRMWIi86SLRRMrGl1dimzbtm2hsoyMDKSnpyM7O1tqJSr4A3rhwgVexpJRWWZjfuWVVxAREQFLS0skJCQYKtQKwegTnn///RetW7eGlZUVfv/9d1y6dAlff/01KleuLNWZN28evv32WyxZsgTHjh2Dg4MDOnXqhKysLBkjJzI/ulg0UZ8q2lD5SpUq4c6dO7hx44bad2K+hw8fSvdzcnIQFhaGefPmcWZlA9J0NuaCPv/8c5w+fbpMl7WM8TNgNMtY6G3qQx0ZP368aNOmTbHbVSqVcHd3F19++aVU9vDhQ2FjYyPWrl2r8XE40zKRZvJnGX5+puH8sudnkDVkXJrMaKspXc5qa8hjFax3+vRpqfz48eMCgHB1dVWbWTk6OlqsXr1a3L17t8zxV2TG9D7R5jNgyJmW9X0ss5lpeevWrXjllVfw9ttvw83NDQEBAVi+fLm0PTExESkpKQgJCZHKnJ2d0aJFCxw5cqTY/T59+hTp6elqNyIqnS4WTdQ1DpUvWp06daT7np6e+Prrr/HJJ5+ozaz85Zdfok+fPvjjjz+ksgcPHmDv3r0aL5RK+nP27Fn079+/1Jm59fEZMJqWGV3ReaqlYzY2NsLGxkZMmDBBnD59WixdulTY2tqKlStXCiGEOHTokAAg7ty5o/Z3b7/9tnjnnXeK3e+0adOKXHOGLTxEmsn/rwqA2LFjh8jNzZUlDk3XkdI2PmP6z12bY2kb07hx40SLFi3ElStXpLJ169YJAKJZs2Zqda9duyays7M1fFYVgz7fJ0+ePBHVq1cXAMTcuXOL3U9ZPgPG9v4uD7Np4VGpVGjSpAlmz56NgIAADBo0CB999BGWLFlSrv1OmDABaWlp0i0pKUlHERNVDMXNMmxoHCpfPl988QWOHj2qNgosKysLvr6+aNGihVQmhEBgYCCcnJxw/vx5tboFZwkm3bG1tcXcuXMBANOnT8e1a9eKrMfPgGaMPuHx8PBAgwYN1Mrq16+PW7duAfj/+Sju3bunVufevXtqc1U8z8bGBk5OTmo3IjI95jBUXpOOpobsjBoeHo4bN27gm2++kcoePHiA7Oxs5ObmqiVHX3/9NapUqYLZs2er7UNweQSdCA8PR3BwMLKystSWcyn4HijLZ8AYOzfrm9EnPK1bty40NO/y5cvw9fUF8GzNGnd3d+zZs0fanp6ejmPHjiEwMNCgsRKR4elrqLyhZraOjo5W+6euc+fO8PPzU+tzoUkdfbCyspLuV61aFampqbh69araul5nzpzBw4cPYWtrK5WlpaXBzc0NHTt2NPsRYfp+nygUCixZsgRWVlb4888/pfKC7wFtPwNyvZ9kp/OLaTp2/PhxYWlpKWbNmiWuXLkiVq9eLezt7cV///tfqc7cuXOFi4uL2LJlizh37pzo1q2b8Pf3F0+ePNH4OBylRaQdQ47yKEl+/4XnR42hhP4LulLe1yB/xFtRMeePeNOkji5j0lZ2drY4c+aMSE5Olsr27t0rAAhfX1+1uhERESIkJETs2rVL73EZE128T4p7bysUCvHrr79q/BnQ5v1kbn14jD7hEUKIbdu2iRdffFHY2NiIevXqiWXLlqltV6lUYsqUKaJ69erCxsZGBAcHi4SEBK2OwYSHSDvGkvAIYbxD5UuiSUdTLy8vvXRG1benT5+KkydPip07d6qVv/TSSwKA2LJli1R29uxZERQUJKZOnWroMA2mPOdE0w7JGzZsKPUzoG3nZiY8ZooJD5F2jOGHtaCi5iDx9vY2ymRHCCH27dtX7A+Ptrd9+/ZJ+zW281LQ6dOnxfLly8X9+/elsh9++EEAEMHBwWp1BwwYIAYMGKD1P6/GqDznRNP3yb59+0r9DGizr/LGravnrwlNf78tQURkBkJDQxESEgJnZ2cAz2a07dixo9EuraDLTtTG3CG7oICAAAQEBKiVderUCatWrZLOG/BsdO769euRmZmJTz75RCr//fff8cMPP6Br167o37+/ocKWlTYdkt97770SPwPm0MG/PJjwEJHZMJah8prQ5Xpjprx2mZeXF/r27atWplKp8MsvvyAuLk5t8sQDBw4gOjoa1apVkxIeIQTefPNN+Pr64vPPP0eVKlUMGb7eadshuaTPQEVfC08hBMcOAs9Gdjk7OyMtLY1D1Ik0kJmZiUqVKgF4toilPkcyacoYYypOXl4e/Pz8kJycXOQQboVCgRo1agBAiXW8vLyQmJgo/bCZ0mugrbi4OOzbtw+NGjVChw4dAAApKSnw8PCAUqlERkaGNILs+++/x86dOxEeHi77avHlOSelvU+AZ/P13Lt3D05OTiUeS5P3XMH3k67eS/p+T2r6+230w9KJiMyRhYUFFi5cCACFFmPNf7xw4cJS6yxYsEDtv3hDDaeXQ+PGjTFq1Cgp2QGePd81a9Zg7ty5asPl9+/fj61btyIxMVEqe/jwIRo1aoS+ffuazLwzJb1P8mVlZeH1119HampqmfdV3PvJrOi895CJYqdlIu0YY+dYY4ypNJp0tja1DtnG4NixY2LhwoXi/PnzUtmBAwcEAOHj46NWd8iQIaJFixZi27ZtUplKpRIqlUonsejifVnce2D27NmicuXKAoBo0KCBSEhIKPVYmr6fzK3TMhOe/2HCQ6QdY0wujDEmTWiyLpmxrF1mylJTU8XWrVvFunXr1MoDAgIEALFp0yap7MSJE8LZ2Vl0795dre7Dhw91lghpq7j3wIULF6QEpuCw85I+A5q8n8wt4eElLSIimWnS2dqUOmQbq8qVK6Nr167o1auXWvmaNWvw66+/onXr1lLZhQsXpLUWC3rjjTfg7OyM3bt3S2VpaWm4du2a3i+TFfceaNiwIQ4dOoS6desiPT29XPvSB2NZxoIJDxERVWj16tXD22+/jWrVqkll7733Hs6dO4d58+ZJZUIIXLlyBY8ePVIbyfT777+jVq1aCAkJUdvvnj17cPbsWWRnZ+v9Ofj6+uLgwYPYtm2bVJabm4vMzEy9H7skxrSMBRMeIioTc+4cS2RjY4OXXnoJr7zyilSmUCiQnJyMixcvom7dulL533//DRsbG7VFVYFnSVPjxo1x4cIFqSwuLg5Lly7F6dOndR5z1apV0bRpU+nxhg0b4O/vj3nz5iEjI0Pr/ZW3ZSY6OhphYWFITk5WK09OTkZYWJjBkx4mPERERBqytrZGgwYN1BZWHTFiBDIzM9VWl3/y5Alq166NKlWqqCVHv/32GwYPHoxvv/1Wbb/Dhw/HrFmz8PDhQ53F+uuvv+Lvv//G+PHj4efnhzlz5uDRo0ca/W15W2by8vIQERFR5PD3/LKRI0ca9PIWEx4iIqJysrCwgKOjo/TYzs4Ohw4dwj///KPWAlqzZk107twZrVq1ksrS09MRGRmJyZMnqw0XX7JkCUJCQrBy5coyxbRhwwasXLkStWrVwoMHDzBx4kT4+flh5syZhfomFaSLlpnY2Fjcvn272O1CCCQlJSE2NlbzJ1ROTHiIiIgMpHfv3ti+fTsGDRoklalUKsyaNQtDhw5VW2Lj+PHj2LNnD5KSkgrtp3Hjxmr9c27evIkbN26otZhYWloiPDwc8fHx+OWXX1C3bl2kpqZiypQpGDhwYJHxadsyk5mZCYVCAYVCoRaPMS5jwYSHiIhIRi4uLpg4cSIiIyPVykeNGoWVK1eiW7duhf7mwYMHai1HM2fOhL+/P+bOnSuVPXnyBEuWLMG+ffvQp08fXLx4EatXr0a9evUwfPhwqV5cXBxGjx6NnTt34o8//tBJy4wxLmPBtbSIiIiM0EsvvYSXXnqpyG2bN29We5ydnQ1ra2vUrFlTKrt+/TqGDBkCFxcXpKamwsLCAr1790Zqaip+++03qd62bdswf/58zJ8/X61vUklKa5kJCgqCl5dXqctYBAUFaXQ8XWALDxERkYkpOBoLAH7++Wc8fvwYPXv2VCvv2rUrOnbsqNY3aNOmTViyZIn0OCgoCD179oSFhQVycnI0On5pLTPGuIwFW3iIyGzkD5UnqogsLCxQuXJltc/A1q1bC9UbNGgQXn75ZSxYsAAA0L59e6hUKmzcuBH16tXD3bt3S+zUrFAo8OGHH+Kll15CcHBwsfVCQ0MRFRWFjz/+WK0DtJeXFxYsWGDwRV25Wvr/cLV0IpKLJqtJm/Mq6KQZXb4Hnt+XSqXChQsXkJOTg3/++QdhYWEAoJY8KRSKEv+hcHJyQmhoKEaOHIlGjRrhyJEjsLKyQvXq1eHj4wMA2LFjBzp27KjTlh1Nf7/ZwkNERFTBOTo6IjAwUHpcXMvMV199hczMTCxatAhKpRLOzs7Yu3cvgGedpFeuXIn+/fsDAIYOHYq4uDi148i5LAoTHiIiIlITGhqKkJAQaZj88y0zH3zwAQD1lqJZs2YhOztbmrCwWrVqcHBwkH15i3xMeIiIiKgQbRcYHTp0qNpltl27diE3NxcnT55Uaz2SCxMeIiKZadLZmh2ySZfvAUO9nywtLYsdWm9oHJZOREREZo8JDxEREZk9JjxERERk9pjwEBERkdljwkNERERmjwkPERERlUleXp50/8CBA2qPjQ0THiIiItJadHS0NMkgAHTu3Bl+fn6Ijo6WMariMeEhIiIirURHRyMsLExt6QkASE5ORlhYmFEmPUx4iIiISGN5eXmIiIgocuLC/LKRI0dKl7fyJzkUQsi66C0THiIiItJYbGwsbt++Xex2IQSSkpIQGxtrwKhKx6UliIiIqJDilp+4e/euRn+vaT1DYQsPERERaczDw0On9QyFCQ8RERFpLCgoCF5eXlAoFEVuVygU8Pb2RlBQkIEjKxkTHiIiItKYhYUFFi5cCACFkp78xwsWLICFhYXBYysJEx4iIiLSSmhoKKKiouDp6alW7uXlhaioKISGhsoUWfEUoqgeSRVQeno6nJ2dkZaWBicnJ7nDISIiMnr5v50AsGPHDnTs2NHgLTua/n6zhYeIiIjKpGBy07ZtW6O7jFUQEx4iIiIye0x4iIiIyOyZXMIzd+5cKBQKjBw5UirLysrCsGHDUKVKFVSqVAk9e/bEvXv35AuSiIiIjIpJJTwnTpzA0qVL8fLLL6uVjxo1Ctu2bcOGDRvw559/4s6dO0bZQ5yIiIjkYTIJT0ZGBvr06YPly5ejcuXKUnlaWhp+/PFHfPPNN+jQoQOaNm2KFStW4PDhwzh69KiMERMREZGxMJmEZ9iwYejSpQtCQkLUyk+dOoWcnBy18nr16sHHxwdHjhwpdn9Pnz5Fenq62o2IiIjMk0ksHrpu3TqcPn0aJ06cKLQtJSUF1tbWcHFxUSuvXr06UlJSit3nnDlzMGPGDF2HSkREREbI6Ft4kpKSEBERgdWrV8PW1lZn+50wYQLS0tKkW1JSks72TURERMbF6BOeU6dO4f79+2jSpAksLS1haWmJP//8E99++y0sLS1RvXp1ZGdn4+HDh2p/d+/ePbi7uxe7XxsbGzg5OandiIiIyDwZ/SWt4OBgnD9/Xq3sgw8+QL169TB+/Hh4e3vDysoKe/bsQc+ePQEACQkJuHXrFgIDA+UImYiIiIyM0Sc8jo6OePHFF9XKHBwcUKVKFal84MCBGD16NFxdXeHk5IQRI0YgMDAQLVu2lCNkIiKiCsHBwQGmsiSn0Sc8mpg/fz6USiV69uyJp0+folOnTvj+++/lDouIiIiMBFdL/x+ulk5ERGR6uFo6ERER0f8w4SEiIiKzx4SHiIiIzB4THiIiIjJ7THiIiIjI7DHhISIiIrPHhIeIiIjMHhMeIiIiMntMeIiIiMjsMeEhIiIis8eEh4iIiMweEx4iIiIye0x4iIiIyOxZyh2AschfND49PV3mSIiIiEhT+b/b+b/jxWHC8z+PHj0CAHh7e8scCREREWnr0aNHcHZ2Lna7QpSWElUQKpUKd+7cgaOjIxQKhdzhGJ309HR4e3sjKSkJTk5OcodD4DkxNjwfxoXnw7jo83wIIfDo0SN4enpCqSy+pw5beP5HqVTCy8tL7jCMnpOTE788jAzPiXHh+TAuPB/GRV/no6SWnXzstExERERmjwkPERERmT0mPKQRGxsbTJs2DTY2NnKHQv/Dc2JceD6MC8+HcTGG88FOy0RERGT22MJDREREZo8JDxEREZk9JjxERERk9pjwEBERkdljwkOSyMhI+Pn5wdbWFi1atMDx48eLrbt8+XIEBQWhcuXKqFy5MkJCQkqsT2WjzTkpaN26dVAoFOjevbt+A6xgtD0fDx8+xLBhw+Dh4QEbGxvUqVMHO3bsMFC05k/b87FgwQLUrVsXdnZ28Pb2xqhRo5CVlWWgaM3bgQMH0LVrV3h6ekKhUGDz5s2l/s3+/fvRpEkT2NjYoFatWli5cqV+gxREQoh169YJa2tr8dNPP4mLFy+Kjz76SLi4uIh79+4VWb93794iMjJSnDlzRsTHx4v+/fsLZ2dncfv2bQNHbr60PSf5EhMTRY0aNURQUJDo1q2bYYKtALQ9H0+fPhWvvPKK6Ny5szh48KBITEwU+/fvF3FxcQaO3Dxpez5Wr14tbGxsxOrVq0ViYqLYuXOn8PDwEKNGjTJw5OZpx44dYtKkSSI6OloAEJs2bSqx/vXr14W9vb0YPXq0uHTpkli0aJGwsLAQMTExeouRCQ8JIYRo3ry5GDZsmPQ4Ly9PeHp6ijlz5mj097m5ucLR0VH8/PPP+gqxwinLOcnNzRWtWrUSP/zwgwgPD2fCo0Pano/FixeLmjVriuzsbEOFWKFoez6GDRsmOnTooFY2evRo0bp1a73GWRFpkvCMGzdONGzYUK2sV69eolOnTnqLi5e0CNnZ2Th16hRCQkKkMqVSiZCQEBw5ckSjfTx+/Bg5OTlwdXXVV5gVSlnPyWeffQY3NzcMHDjQEGFWGGU5H1u3bkVgYCCGDRuG6tWr48UXX8Ts2bORl5dnqLDNVlnOR6tWrXDq1Cnpstf169exY8cOdO7c2SAxk7ojR46onT8A6NSpk8a/OWXBxUMJ//zzD/Ly8lC9enW18urVq+Ovv/7SaB/jx4+Hp6dnoTcwlU1ZzsnBgwfx448/Ii4uzgARVixlOR/Xr1/H3r170adPH+zYsQNXr17F0KFDkZOTg2nTphkibLNVlvPRu3dv/PPPP2jTpg2EEMjNzcXgwYMxceJEQ4RMz0lJSSny/KWnp+PJkyews7PT+THZwkPlNnfuXKxbtw6bNm2Cra2t3OFUSI8ePULfvn2xfPlyVK1aVe5wCIBKpYKbmxuWLVuGpk2bolevXpg0aRKWLFkid2gV0v79+zF79mx8//33OH36NKKjo7F9+3Z8/vnncodGBsIWHkLVqlVhYWGBe/fuqZXfu3cP7u7uJf7tV199hblz5+KPP/7Ayy+/rM8wKxRtz8m1a9dw48YNdO3aVSpTqVQAAEtLSyQkJOCFF17Qb9BmrCyfEQ8PD1hZWcHCwkIqq1+/PlJSUpCdnQ1ra2u9xmzOynI+pkyZgr59++LDDz8EALz00kvIzMzEoEGDMGnSJCiV/P/fkNzd3Ys8f05OTnpp3QHYwkMArK2t0bRpU+zZs0cqU6lU2LNnDwIDA4v9u3nz5uHzzz9HTEwMXnnlFUOEWmFoe07q1auH8+fPIy4uTrq99dZbePXVVxEXFwdvb29Dhm92yvIZad26Na5evSolngBw+fJleHh4MNkpp7Kcj8ePHxdKavKTUcElJQ0uMDBQ7fwBwO7du0v8zSk3vXWHJpOybt06YWNjI1auXCkuXbokBg0aJFxcXERKSooQQoi+ffuKTz/9VKo/d+5cYW1tLaKiosTdu3el26NHj+R6CmZH23PyPI7S0i1tz8etW7eEo6OjGD58uEhISBC//fabcHNzEzNnzpTrKZgVbc/HtGnThKOjo1i7dq24fv262LVrl3jhhRfEO++8I9dTMCuPHj0SZ86cEWfOnBEAxDfffCPOnDkjbt68KYQQ4tNPPxV9+/aV6ucPSx87dqyIj48XkZGRHJZOhrNo0SLh4+MjrK2tRfPmzcXRo0elbe3atRPh4eHSY19fXwGg0G3atGmGD9yMaXNOnseER/e0PR+HDx8WLVq0EDY2NqJmzZpi1qxZIjc318BRmy9tzkdOTo6YPn26eOGFF4Stra3w9vYWQ4cOFf/++6/hAzdD+/btK/I3If8chIeHi3bt2hX6m8aNGwtra2tRs2ZNsWLFCr3GqBCCbXlERERk3tiHh4iIiMweEx4iIiIye0x4iIiIyOwx4SEiIiKzx4SHiIiIzB4THiIiIjJ7THiIiIjI7DHhISIiIrPHhIeIiIjMHhMeIiIiMntMeIiIiMjsMeEhIrO0du1a2NnZ4e7du1LZBx98gJdffhlpaWkyRkZEcuDioURkloQQaNy4Mdq2bYtFixZh2rRp+Omnn3D06FHUqFFD7vCIyMAs5Q6AiEgfFAoFZs2ahbCwMLi7u2PRokWIjY1lskNUQbGFh4jMWpMmTXDx4kXs2rUL7dq1kzscIpIJ+/AQkdmKiYnBX3/9hby8PFSvXl3ucIhIRmzhISKzdPr0abRv3x5Lly7FypUr4eTkhA0bNsgdFhHJhH14iMjs3LhxA126dMHEiRPx3nvvoWbNmggMDMTp06fRpEkTucMjIhmwhYeIzEpqaipatWqF9u3bY8mSJVJ5ly5dkJeXh5iYGBmjIyK5MOEhIiIis8dOy0RERGT2mPAQERGR2WPCQ0RERGaPCQ8RERGZPSY8REREZPaY8BAREZHZY8JDREREZo8JDxEREZk9JjxERERk9pjwEBERkdljwkNERERm7/8AqtSV7Vtjp6gAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def poisson_interval(k, alpha=0.32): \n",
" \"\"\" Uses chi2 to get the poisson interval.\n",
" \"\"\"\n",
" a = alpha\n",
" low, high = (chi2.ppf(a/2, 2*k) / 2, chi2.ppf(1-a/2, 2*k + 2) / 2)\n",
" if k == 0: \n",
" low = 0.0\n",
" return k - low, high - k\n",
"\n",
"x = np.linspace(0.1, 1, 50) # Dependent variable\n",
"\n",
"# Mean expected counts\n",
"y_mu = bump_forward_model(x, \n",
" amp_s=50, mu_s=0.8, std_s=0.05, # Signal params\n",
" amp_b=50, exp_b=-0.5) # Background params\n",
"\n",
"# Mean expected counts without signal\n",
"y_mu_sig = bump_forward_model(x, \n",
" amp_s=0., mu_s=0.8, std_s=0.05, # Signal params\n",
" amp_b=50, exp_b=-0.5) # Background params\n",
"\n",
"\n",
"# Realized counts\n",
"y = onp.random.poisson(y_mu)\n",
"y_err = np.array([poisson_interval(k) for k in y.T]).T\n",
"\n",
"# Plot\n",
"plt.plot(x, y_mu, color='k', ls='--', label=\"Mean expected counts\")\n",
"plt.plot(x, y_mu_sig, color='k', ls=':', label=\"Mean expected counts, bkg\")\n",
"plt.errorbar(x, y, yerr=y_err, fmt='o', color='k', label=\"Realized counts\")\n",
"\n",
"plt.xlabel(\"$x$\")\n",
"plt.ylabel(\"$y$\")\n",
"\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2f2f049d-1244-42f4-917b-f721afd9c98c",
"metadata": {},
"outputs": [],
"source": [
"import numpyro\n",
"import numpyro.distributions as dist\n",
"from numpyro.infer.autoguide import AutoMultivariateNormal\n",
"\n",
"from tinygp import kernels, GaussianProcess\n",
"\n",
"jax.config.update(\"jax_enable_x64\", True)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ba75fd8a-75d0-4316-8782-f540fa732ae1",
"metadata": {},
"outputs": [],
"source": [
"def model(x, y):\n",
" \n",
" # The parameters of the GP model\n",
" mean = numpyro.param(\"mean\", np.zeros(()))\n",
" sigma = numpyro.param(\"sigma\", np.ones(()), constraint=dist.constraints.positive)\n",
" rho = numpyro.param(\"rho\", np.ones(()), constraint=dist.constraints.positive)\n",
" \n",
" amp_b = numpyro.sample(\"amp_b\", dist.Uniform(0., 200.))\n",
"\n",
" # Set up the kernel and GP objects\n",
" kernel = sigma ** 2. * kernels.ExpSquared(rho)\n",
" gp = GaussianProcess(kernel, x, diag=1e-5, mean=np.zeros(()))\n",
"\n",
" # This parameter has shape (num_data,) and it encodes our beliefs about\n",
" # the process rate in each bin\n",
" log_rate = numpyro.sample(\"log_rate\", gp.numpyro_dist())\n",
" \n",
" # Mean expected counts without signal\n",
" rate_bkg = bump_forward_model(x, \n",
" amp_s=0., mu_s=0.8, std_s=0.05, # Signal params\n",
" amp_b=amp_b, exp_b=-0.5) # Background params\n",
" \n",
" rate_tot = np.exp(log_rate) + rate_bkg\n",
" \n",
" numpyro.sample(\"obs\", dist.Poisson(rate_tot), obs=y)\n",
"\n",
"guide = AutoMultivariateNormal(model)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "20d109a5-3acd-4dc6-a9ca-488150002a9a",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█| 30000/30000 [00:39<00:00, 756.54it/s, init loss: 2162058.5566, avg. loss\n"
]
}
],
"source": [
"optim = numpyro.optim.Adam(1e-1)\n",
"svi = numpyro.infer.SVI(model, guide, optim, numpyro.infer.Trace_ELBO(16))\n",
"results = svi.run(jax.random.PRNGKey(42), 30000, x, y, progress_bar=True)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5bdb0823-48aa-49ca-b801-aab581e86bb9",
"metadata": {},
"outputs": [],
"source": [
"num_samples = 10000\n",
"\n",
"samples = guide.sample_posterior(\n",
" rng_key=jax.random.PRNGKey(4),\n",
" params=results.params,\n",
" sample_shape=(num_samples,)\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1e41196c-104a-4bf0-83a2-898ac4e0cb51",
"metadata": {},
"outputs": [],
"source": [
"# Mean expected counts without signal\n",
"samples_bkg = np.array([bump_forward_model(x, \n",
" amp_s=0., mu_s=0.8, std_s=0.05, # Signal params\n",
" amp_b=amp_b, exp_b=-0.5) # Background params\n",
" for amp_b in samples['amp_b']])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "dba5ff7c-3631-48d6-b3f5-607cf2f4088e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x16b6d62b0>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGwCAYAAACtlb+kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC83UlEQVR4nOzdd3hUVfrA8e+9UzMlvZOE3pEOCgqCoIJl1ciu7ror/nTtDV1XZK3YsCuWxV3Xvrqssqi7iCjYQEAEJEgNLZCQXieZ3u7vjyGTDGkTSOd8nmceMvfeufdMmMy8c8573iMpiqIgCIIgCILQg8md3QBBEARBEIT2JgIeQRAEQRB6PBHwCIIgCILQ44mARxAEQRCEHk8EPIIgCIIg9Hgi4BEEQRAEoccTAY8gCIIgCD2eurMb0FX4/X4KCgowm81IktTZzREEQRAEIQyKolBTU0Nqaiqy3HQ/jgh4jikoKCA9Pb2zmyEIgiAIwgnIy8sjLS2tyf0i4DnGbDYDgV9YZGRkJ7dGEARBEIRwVFdXk56eHvwcb4oIeI6pHcaKjIwUAY8gCIIgdDMtpaOIpGVBEARBEHq8Tg941q5dy8UXX0xqaiqSJPHpp5+G7JckqdHbs88+GzymT58+DfY/9dRTHfxMBEEQBEHoqjo94LHZbIwaNYrXXnut0f2FhYUht7feegtJkrj88stDjnv00UdDjrv99ts7ovmCIAiCIHQDnZ7DM3v2bGbPnt3k/uTk5JD7n332GdOnT6dfv34h281mc4NjBUEQhLbn8/nweDyd3QzhFKHRaFCpVCd9nk4PeFqjuLiYzz//nHfffbfBvqeeeorHHnuMjIwMfve733HXXXehVjf99FwuFy6XK3i/urq6XdosCILQUyiKQlFREVVVVZ3dFOEUEx0dTXJy8knVyetWAc+7776L2WwmMzMzZPsdd9zB2LFjiY2NZcOGDSxYsIDCwkJeeOGFJs+1aNEiFi5c2N5NFgRB6DFqg53ExEQMBoMo0iq0O0VRsNvtlJSUAJCSknLC55IURVHaqmEnS5IkPvnkEy699NJG9w8ZMoRzzz2XV155pdnzvPXWW9x4441YrVZ0Ol2jxzTWw5Oeno7FYhHT0gVBEI7j8/nYt28fiYmJxMXFdXZzhFNMeXk5JSUlDBo0qMHwVnV1NVFRUS1+fnebHp5169aRnZ3Nv//97xaPPf300/F6vRw+fJjBgwc3eoxOp2syGBIEQRBC1ebsGAyGTm6JcCqqfd15PJ4Tzufp9Fla4XrzzTcZN24co0aNavHYrKwsZFkmMTGxA1omCIJw6hDDWEJnaIvXXaf38FitVg4cOBC8n5OTQ1ZWFrGxsWRkZACB7qqPP/6Y559/vsHjN27cyKZNm5g+fTpms5mNGzdy11138fvf/56YmJgOex6CIAiCIHRdnR7wbNmyhenTpwfv33333QDMnTuXd955B4ClS5eiKAq//e1vGzxep9OxdOlSHnnkEVwuF3379uWuu+4KnkcQBEEQBKFLJS13pnCTngRBEE5FTqeTnJwc+vbti16v7+zmdDl9+vRh3rx5zJs3r7Ob0mbeeecd5s2b1yXKEDT3+gv387vb5PAIgiAIwom45pprQpYeiouLY9asWfzyyy+d3TShA4mARxAEQejxZs2aFVx66Ouvv0atVnPRRRd1drOa5Xa7O7sJPYoIeNqRzWYLfqOw2Wyd3RxBEIS2pSjgtnXOrZXZGDqdjuTkZJKTkxk9ejT33XcfeXl5lJaWAjB//nwGDRqEwWCgX79+PPjggw2Wz/jf//7HhAkT0Ov1xMfHc9lllzV5vX/84x9ER0fz9ddfA1BTU8NVV12F0WgkJSWFF198kWnTpoUMgfXp04fHHnuMq6++msjISG644QYA/vOf/zB8+HB0Oh19+vRpMIGnsYW3o6Ojg3mwhw8fRpIkli9fzvTp0zEYDIwaNYqNGzeGPOadd94hIyMDg8HAZZddRnl5edi/3+6g05OWBUEQhG7KY4cnUzvn2n8pAK3xhB5qtVr55z//yYABA4JFFM1mM++88w6pqans2LGD66+/HrPZzL333gvA559/zmWXXcb999/Pe++9h9vtZuXKlY2e/5lnnuGZZ57hq6++YuLEiUBgQs769ev573//S1JSEg899BA///wzo0ePDnnsc889x0MPPcTDDz8MwNatW/nNb37DI488whVXXMGGDRu45ZZbiIuL45prrmnV877//vt57rnnGDhwIPfffz+//e1vOXDgAGq1mk2bNnHdddexaNEiLr30UlatWhVsQ08hAh5BEAShx1uxYgUmkwkI9L6npKSwYsUKZDkw0PHAAw8Ej+3Tpw/33HMPS5cuDQY8TzzxBFdeeWXIkkSN1YWbP38+77//Pt9//z3Dhw8HAr077777Lh9++CEzZswA4O233yY1tWGweM455/CnP/0peP+qq65ixowZPPjggwAMGjSI3bt38+yzz7Y64Lnnnnu48MILAVi4cCHDhw/nwIEDDBkyhMWLFzNr1qzg8x00aBAbNmxg1apVrbpGVyYCHkEQBOHEaAyBnpbOunYrTJ8+nSVLlgBQWVnJX//6V2bPns1PP/1E7969+fe//83LL7/MwYMHsVqteL3ekBk/WVlZXH/99c1e4/nnn8dms7Flyxb69esX3H7o0CE8Hk+wtwcgKiqq0ZUAxo8fH3J/z549XHLJJSHbzjzzTF566SV8Pl+rqg6PHDky+HPtmlQlJSUMGTKEPXv2NBiimzRpUo8KeEQOjyAIgnBiJCkwrNQZt1ZW3jUajQwYMIABAwYwYcIE/vGPf2Cz2XjjjTfYuHEjV111FRdccAErVqxg27Zt3H///SFJwxERES1eY8qUKfh8Pj766KNW/yrrt7O1JEni+Aozx+cfAWg0mpDHAPj9/lZfr7sSAY8gCIJwypEkCVmWcTgcbNiwgd69e3P//fczfvx4Bg4cyJEjR0KOHzlyZDABuSkTJ07kiy++4Mknn+S5554Lbu/Xrx8ajYbNmzcHt1ksFvbt29diO4cOHcr69etDtq1fvz5kEc2EhAQKCwuD+/fv34/dbm/x3MdfZ9OmTSHbfvzxx1ado6sTQ1qCIAhCj+dyuSgqKgICQ1qvvvoqVquViy++mOrqanJzc1m6dCkTJkzg888/55NPPgl5/MMPP8yMGTPo378/V155JV6vl5UrVzJ//vyQ4yZPnszKlSuZPXs2arWaefPmYTabmTt3Ln/+85+JjY0lMTGRhx9+GFmWW1wj6k9/+hMTJkzgscce44orrmDjxo28+uqr/PWvfw0ec8455/Dqq68yadIkfD4f8+fPD+nNCccdd9zBmWeeyXPPPccll1zCl19+2aOGs0D08AiCIAingFWrVpGSkkJKSgqnn346mzdv5uOPP2batGn86le/4q677uK2225j9OjRbNiwIZgkXGvatGl8/PHH/Pe//2X06NGcc845/PTTT41e66yzzuLzzz/ngQce4JVXXgHghRdeYNKkSVx00UXMnDmTM888k6FDh7ZYtXrs2LF89NFHLF26lBEjRvDQQw/x6KOPhiQsP//886SnpzNlyhR+97vfcc8997R6VfszzjiDN954g8WLFzNq1Ci++uqrkETunkAsLXFMeywtYbPZgrMCrFbrCY3NCoIgdAViaYm2ZbPZ6NWrF88//zzXXXddZzeny2uLpSXEkJYgCIIgtLNt27axd+9eJk6ciMVi4dFHHwVoMANLaD8i4BEEQRCEDvDcc8+RnZ2NVqtl3LhxrFu3jvj4+M5u1ilDBDwdxOf3dXYTBEEQhE4yZswYtm7d2tnNOKWJpOUOUuYo6+wmCIIgCMIpSwQ8HaTc0bMWYRMEQRCE7kQEPB3E4rZg97SuEJQgCIIgCG1DBDwdxOa0YHFZOrsZgiAIgnBKEgFPO3L5XMGfLdYCShwlDdY7EQRBEASh/YmApx3pVLrgz/tKd2FxVmH1WDuxRYIgCIJwahIBTwfZVpWN2yWGtQRBEDpLUVERd955JwMGDECv15OUlMSZZ57JkiVLgott9unTB0mSkCQJo9HI2LFj+fjjjzu55UJbEAFPBznsLMFqK6HEXoJf8Xd2cwRBEE4phw4dYsyYMXz11Vc8+eSTbNu2jY0bN3LvvfeyYsUK1qxZEzz20UcfpbCwkG3btjFhwgSuuOIKNmzY0ImtF9qCKDzYnjzOkLtZxVuJie5PjbuGKF1UJzVKEATh1HPLLbegVqvZsmVLyLqG/fr145JLLgnJrzSbzSQnJ5OcnMxrr73GP//5T/73v/8xefLkzmi60EZEwNOeNKELnG0s38VZLgtVzioR8AiC0O0pioLD0zlV5CM0KiRJCuvY8vLyYM9OU4s4N3UutVqNRqPB7XafcFuFrkEEPO3I56t7I3DttZJ/mkRZTR5Rkb3oZe6FWha/fkEQui+Hx8ewh77slGvvfvR8DNrw3kMPHDiAoigMHjw4ZHt8fDxOZ6An/tZbb+Xpp58O2e92u3n++eexWCycc845bdNwodOIHJ52snz5coYNGxa8f/DFXLL/lM3S5Z9gdVqodld3YusEQRCEn376iaysLIYPH47LVVdGZP78+ZhMJgwGA08//TRPPfUUF154YSe2VGgLoouhHSxfvpw5c+Y0qLnjrfTy2ZPfMjp5Ahm/70usPraTWigIgnDyIjQqdj96fqddO1wDBgxAkiSys7NDtvfr1y9wroiIkO1//vOfueaaazCZTCQlJYU9dCZ0bSLgaWM+n48777yz2QKDrzz+N8679BIyIjPQqDQd2DpBEIS2I0lS2MNKnSkuLo5zzz2XV199ldtvv73JPJ5a8fHxDBgwoINaJ3QUMaTVxtatW8fRo0ebPaaixMKPP/yAxS1q8giCIHSEv/71r3i9XsaPH8+///1v9uzZQ3Z2Nv/85z/Zu3cvKlX4PUZC99T1Q/NuprCwMKzjSvPzKXeUEx8R384tEgRBEPr378+2bdt48sknWbBgAUePHkWn0zFs2DDuuecebrnlls5uotDORMDTxlJSUsI6zqGtotxRjtPrRK/Wt/wAQRAE4aSkpKTwyiuv8MorrzR5zOHDhzuuQUKHEkNabWzKlCmkpaU1m+SmidVQkVqF3VkllpoQBEEQhA4gAp42plKpWLx4MdCwkJV07Jb8u2R+sR/G5aygzFHW8Y0UBEEQhFOMCHjaQWZmJsuWLSM1NTVke2pSLB/9xsDE0yLwKX72FG+nsKIguFCdzWbrpBYLgiAIQs8mAp52kpmZye7du4P3V65cycHNa5hx9gguPBbYbKrYgVMMaQmCIAhCuxMBTzuqP81x6tSp6GJ64c2YxGyrHUlR2GvLp8ZW0oktFARBEIRTgwh4OlJEDNr+5xAvaRjvDJQx31WyvZMbJQiCIAg9nwh4OpJKgzl+CJaUEVxoDQxr/VSxu4UHCYIgCIJwskTA08FkYzxKnymca7ejURTy3RXBfX7F34ktEwRBEISeq9MDnrVr13LxxReTmpqKJEl8+umnIfuvueaa4Cym2tusWbNCjqmoqOCqq64iMjKS6OhorrvuOqxWawc+i1bQRxPRZypabSRn2x0hu2rcNZ3UKEEQBEHo2To94LHZbIwaNYrXXnutyWNmzZpFYWFh8Pavf/0rZP9VV13Frl27WL16NStWrGDt2rXccMMN7d30EyPLGOMGUNNrTHBYq1aFs6KJBwmCIAjdVWNf5us7fPgwkiSRlZXV6W3pyTp9aYnZs2cze/bsZo/R6XQkJyc3um/Pnj2sWrWKzZs3M378eABeeeUVLrjgAp577rkGtXBquVwuXC5X8H51dfUJPoMTEBGL1P8cphz6DlNk3TBWka2Ifp5+GDSGEzqtzWbDZDIBYLVaW1wRWBAE4VRSVFTEE088weeff05+fj6JiYmMHj2aefPmMWPGjE5rV3p6OoWFhcTHi7UV21On9/CE47vvviMxMZHBgwdz8803U15eHty3ceNGoqOjg8EOwMyZM5FlmU2bNjV5zkWLFhEVFRW8paent+tzCKEzYU4/A485hXNs9uBmh8tKuaO8mQcKgiAIJ+Lw4cOMGzeOb775hmeffZYdO3awatUqpk+fzq233tqpbVOpVCQnJ6NWd3ofRI/W5QOeWbNm8d577/H111/z9NNP8/333zN79mx8Ph8QiNgTExNDHqNWq4mNjaWoqKjJ8y5YsACLxRK85eXltevzOJ4uKgN72gRm1Qt4tF4PBdYCXD5Xg+NtNpuoyCwIgnCCbrnlFiRJ4qeffuLyyy9n0KBBDB8+nLvvvpsff/wRgNzcXC655BJMJhORkZH85je/obi4OHiORx55hNGjR/PWW2+RkZGByWTilltuwefz8cwzz5CcnExiYiJPPPFEg+sXFhYye/ZsIiIi6NevH8uWLQvuO35I67vvvkOSJL7++mvGjx+PwWBg8uTJZGdnh5zzs88+Y+zYsej1evr168fChQvxer3B/fv372fq1Kno9XqGDRvG6tWr2/JX2u10+XDyyiuvDP582mmnMXLkSPr378933313Ul2QOp0OnU7XFk08MYZY1APOZUzWf4ObDpRtZ0h0DOWOclJNjQ/FCYIgdBWKouDwOlo+sB1EqCOaXaS5voqKClatWsUTTzzR6FB/dHQ0fr8/GOx8//33eL1ebr31Vq644gq+++674LEHDx7kiy++YNWqVRw8eJA5c+Zw6NAhBg0axPfff8+GDRu49tprmTlzJqeffnrwcQ8++CBPPfUUixcv5v333+fKK69kx44dDB06tMl233///Tz//PMkJCRw0003ce2117J+/XoA1q1bx9VXX83LL7/MlClTOHjwYDB39eGHH8bv95OZmUlSUhKbNm3CYrEwb968sH5fPVWXD3iO169fP+Lj4zlw4AAzZswgOTmZkpLQasVer5eKioom8366BLUOc+oYSuL6AoHig6sKNjCq71kUWAtIMCSgkTWd20ZBEIRmOLwOTv/w9JYPbAebfrcp7HzHAwcOoCgKQ4YMafKYr7/+mh07dpCTkxNMcXjvvfcYPnw4mzdvZsKECQD4/X7eeustzGYzw4YNY/r06WRnZ7Ny5UpkWWbw4ME8/fTTfPvttyEBz69//Wv++Mc/AvDYY4+xevVqXnnlFf7617822aYnnniCs88+G4D77ruPCy+8EKfTiV6vZ+HChdx3333MnTsXCHw2PvbYY9x77708/PDDrFmzhr179/Lll18Gc1mffPLJFnNme7IuP6R1vKNHj1JeXk5KSgoAkyZNoqqqiq1btwaP+eabb/D7/SEvtq5IY07BnT4xeP+Qs5ijFXuxuCxUOMSMLUEQhLagKEqLx+zZs4f09PSQfM5hw4YRHR3Nnj17gtv69OmD2WwO3k9KSmLYsGHIshyy7fgv4pMmTWpwv/55GzNy5Mjgz7WfebXn3b59O48++igmkyl4u/766yksLMRutwefT/2JO8e34VTT6T08VquVAwcOBO/n5OSQlZVFbGwssbGxLFy4kMsvv5zk5GQOHjzIvffey4ABAzj//PMBGDp0KLNmzeL666/n9ddfx+PxcNttt3HllVc2OUOry9BHo+s/A3gjuOm/uau5KX4YhbZC4iPiUcmqph8vCILQiSLUEWz6XdOTQ9r72uEaOHAgkiSxd+/ek76uRhPa8y5JUqPb/P6TLyRb/7y1w3e157VarSxcuJDMzMwGj9Pr9Sd97Z6o0wOeLVu2MH369OD9u+++G4C5c+eyZMkSfvnlF959912qqqpITU3lvPPO47HHHgvJv/nggw+47bbbmDFjBrIsc/nll/Pyyy93+HNpNZUac3zd+K2swA5rLmWVOfhUGipdlcRHiGmKgiB0TZIknXAZjY4UGxvL+eefz2uvvcYdd9zRII+nqqqKoUOHkpeXR15eXrCXZ/fu3VRVVTFs2LCTbsOPP/7I1VdfHXJ/zJgxJ3y+sWPHkp2dzYABAxrdX/t8CgsLg71DtcnZp6pOD3imTZvWbHfjl19+2eI5YmNj+fDDD9uyWR1GNsYFf55ttfGF3syK3C/5Q0w/imxFxOnjwk7MEwRBEBr32muvceaZZzJx4kQeffRRRo4cidfrZfXq1SxZsoTdu3dz2mmncdVVV/HSSy/h9Xq55ZZbOPvss0PKnpyojz/+mPHjx3PWWWfxwQcf8NNPP/Hmm2+e8PkeeughLrroIjIyMpgzZw6yLLN9+3Z27tzJ448/zsyZMxk0aBBz587l2Wefpbq6mvvvv/+kn0d31u1yeHocfVTwx6urq5GAzZb9WK0FlDvKsbgsndc2QRCEHqJfv378/PPPTJ8+nT/96U+MGDGCc889l6+//polS5YgSRKfffYZMTExTJ06lZkzZ9KvXz/+/e9/t8n1Fy5cyNKlSxk5ciTvvfce//rXv06q5+j8889nxYoVfPXVV0yYMIEzzjiDF198kd69ewMgyzKffPIJDoeDiRMn8sc//rHR6fKnEkkJJ5vrFFBdXU1UVBQWi4XIyMg2OWc4lY9Djllg5oH0JNYY9JwVO4Jfj7iO1Kg0hsQOaf25RKVlQRDakNPpJCcnh759+4ocEaHDNff6C/fzW/TwdCFWUwLXV1YCsL5iFy57GSX2EqrdHbjshSAIgiD0QCLgaUdGoxFFUVAUpcneltqK0QCr7P0Z7HQz2elFQeGrI6vw+DyU2ko7qsmCIAiC0COJgKcTLV++PGQM95onVtP7ZTsD1hYA8F3ZNrxOC0X2Iuwee1OnEQRBEAShBSLg6STLly9nzpw55Ofnh2wvsPiY/2E1cZsseBU/3x5ZjcPrEIuKCoIgCMJJEAFPJ/D5fNx5552NTsev3bLnw3wUv8Ka4k3gdlBoK+zYRgqCIAhCDyICnk6wbt06jh492uR+BSix+IndY8fl9/BD7jdY3daOa6AgCIIg9DAi4OkEhYXh9daceTQwjLWq8Ackv7c9myQIgiAIPZoIeDpBbZnvlszUeslQVNh8Trbk/9DOrRIEQRCEnksEPJ1gypQppKWlNblkhCRJpCZEMi1DxfXlgZVxvzy6Nrjf4/d0SDsFQRDag81mQ5IkJEnCZrN1dnOEU4QIeDqBSqVi8eLFAA2Cntr7zzx4O9VJg7iwxkYSaqp9ddPSy+2Nz9iqX9Nn7dq1IfcFQRCE9nXNNddw6aWXBu9PmzaNefPmtes1JUni008/bddr9BQi4OkkmZmZLFu2jNTU1JDtaWlpLFu2jKt+dy22ftPQAP9XURFyzFHrUdw+d8i242v6XHDBBfTp04fly5e323MQBEE4ER395eyaa64J9ihpNBr69u3Lvffei9PpbNfrLl++nMcee6xdr9EdHD58GEmSyMrK6tR2iICnE2VmZrJ79+7g/ZUrV5KTk0NmZiaYEjH3noI1MoXLqy1E+uv+q9b/sJ6C6oLg/aZq+uTn5zNnzhwR9AiC0GV01pezWbNmUVhYyKFDh3jxxRf529/+xsMPP9yu14yNjcVsNrfrNYTwiYCnk6lUquDPU6dOrbuvMxEZNwhLv7NZudvN3vv2BI+7/5r7OX346fzr4381X9Pn2LZ58+aJ4S1BEDpdR3w58/l8bNmyhS1btoS87+l0OpKTk0lPT+fSSy9l5syZrF69Orjf7/ezaNEi+vbtS0REBKNGjWLZsmUh573uuuuC+wcPHhxMTWhK/SGt7777LtjLVP92zTXXBI//7LPPGDt2LHq9nn79+rFw4UK83roZuvv372fq1Kno9XqGDRsW0v6m+P1+nnnmGQYMGIBOpyMjIyNk1fQdO3ZwzjnnEBERQVxcHDfccANWa10ZlMaG5S699NKQdvfp04cnn3ySa6+9FrPZTEZGBn//+9+D+/v27QvAmDFjkCSJadOmBX8nEydOxGg0Eh0dzZlnnsmRI0dafE4nSgQ8XZgUlcqaw2bmfOSgyhIasJQUlXDVFVfxxBNPNF/TR1HIy8tj3bp17d1cQRCEJnWlL2c7d+5kw4YNaLXa4LZFixbx3nvv8frrr7Nr1y7uuusufv/73/P9998DgcAhLS2Njz/+mN27d/PQQw/xl7/8hY8++iisa06ePJnCwsLg7ZtvvkGv1zN16lQgUJ/t6quv5s4772T37t387W9/45133gkGJ36/n8zMTLRaLZs2beL1119n/vz5LV53wYIFPPXUUzz44IPs3r2bDz/8kKSkJCCQPH7++ecTExPD5s2b+fjjj1mzZg233XZbq36fAM8//zzjx49n27Zt3HLLLdx8881kZ2cD8NNPPwGwZs0aCgsLWb58OV6vl0svvZSzzz6bX375hY0bN3LDDTc0OZmnTSiCoiiKYrFYFECxWCwdel2r1aoQqDWoWK3WkH1er1dJS0kK7m9wk1BiYmOa3l/v9uGHH3bo8xIEoWdxOBzK7t27FYfDcUKP//bbb8N6r/r2229Pqp1er1fZvHmzsnnzZsXr9SqKoihz585VVCqVYjQaFZ1OpwCKLMvKsmXLFEVRFKfTqRgMBmXDhg0h57ruuuuU3/72t01e69Zbb1Uuv/zy4P25c+cql1xySfD+2Wefrdx5550NHldWVqb069dPueWWW4LbZsyYoTz55JMhx73//vtKSkqKoiiK8uWXXypqtVrJz88P7v/iiy8UQPnkk08abV91dbWi0+mUN954o9H9f//735WYmJiQz57PP/9ckWVZKSoqavI5XHLJJcrcuXOD93v37q38/ve/D973+/1KYmKismTJEkVRFCUnJ0cBlG3btgWPKS8vVwDlu+++a7Rtx2vu9Rfu57e6/UIp4WStW7eOo4XFTR+gQGVFZVjnCrf2jyAIQnsIt+BquMe11vTp01myZAk2m40XX3wRtVrN5ZdfDsCBAwew2+2ce+65IY9xu92MGTMmeP+1117jrbfeIjc3F4fDgdvtZvTo0a1qh8fj4fLLL6d3794hQ2Lbt29n/fr1IcNNPp8Pp9OJ3W5nz549pKenh0x0mTRpUrPX2rNnDy6XixkzZjS5f9SoURiNxuC2M888E7/fT3Z2drAnKBwjR44M/ixJEsnJyZSUlDR5fGxsLNdccw3nn38+5557LjNnzuQ3v/lNu35WiYCnCwv3Dz8mJoaqqqpGu4olSSItLY0pU6a0dfMEQRDCFu4HWXt94BmNRgYMGADAW2+9xahRo3jzzTe57rrrgjkrn3/+Ob169Qp5nE6nA2Dp0qXcc889PP/880yaNAmz2cyzzz7Lpk2bWtWOm2++mby8PH766SfU6rqPYKvVysKFCwOTVo6j1+tbdY1aERERJ/S4+mRZbvDZ4vE0rAWn0WhC7kuShN/vb/bcb7/9NnfccQerVq3i3//+Nw888ACrV6/mjDPOOOl2N0bk8HRh4f7hX3XDVUDTNX1eeumlkORoQRCEjhZOwdX09PQO+XImyzJ/+ctfeOCBB3A4HAwbNgydTkdubi4DBgwIuaWnpwOwfv16Jk+ezC233MKYMWMYMGAABw8ebNV1X3jhBT766CM+++wz4uLiQvaNHTuW7OzsBtcfMGAAsiwzdOhQ8vLyQr4I//jjj81eb+DAgURERPD11183un/o0KFs3749pPjj+vXrkWWZwYMHA5CQkBByTZ/Px86dO1v1vGtzpRrLzxozZgwLFixgw4YNjBgxgg8//LBV524NEfB0YS29QQDEJ8dx6U2X8vaHbzdZ06exbwyCIAgdKZyCqx355ezXv/41KpWK1157DbPZzD333MNdd93Fu+++y8GDB/n555955ZVXePfdd4FA8LBlyxa+/PJL9u3bx4MPPsjmzZvDvt6aNWu49957efbZZ4mPj6eoqIiioiIsFgsADz30EO+99x4LFy5k165d7Nmzh6VLl/LAAw8AMHPmTAYNGsTcuXPZvn0769at4/7772/2mnq9nvnz53Pvvffy3nvvcfDgQX788UfefPNNAK666ir0ej1z585l586dfPvtt9x+++384Q9/CA5nnXPOOXz++ed8/vnn7N27l5tvvpmqqqpW/a4TExOJiIhg1apVFBcXY7FYyMnJYcGCBWzcuJEjR47w1VdfsX//foYOHdqqc7eGCHi6sObeIGqdds0QZJXM2Blj2bmrLuoOqekjCILQBbRUcLUj36/UajW33XYbzzzzDDabjccee4wHH3yQRYsWMXToUGbNmsXnn38enFJ94403kpmZyRVXXMHpp59OeXk5t9xyS9jX++GHH/D5fNx0002kpKQEb3feeScA559/PitWrOCrr75iwoQJnHHGGbz44ov07t0bCPRKffLJJzgcDiZOnMgf//jHkHyfpjz44IP86U9/4qGHHmLo0KFcccUVwdwag8HAl19+SUVFBRMmTGDOnDnMmDGDV199Nfj4a6+9lrlz53L11Vdz9tln069fP6ZPnx7284bA7/rll1/mb3/7G6mpqVxyySUYDAb27t3L5ZdfzqBBg7jhhhu49dZbufHGG1t17taQlMYSP05B1dXVREVFYbFYiIyM7LDr2mw2TCYTEBjDrZ88Vmv58uXccfvt5BfUFRtMiZLR/iEN8/hIHhr3Z+Ii0+gX0Y++iX2bPZcgCMKJcDqd5OTk0Ldv3xPOKalV+34LgS9n5513Xpv17Ph8PrZt2wYEhkvEcH7P0NzrL9zPb9HD0w1kZmaye09d4cGPb+hD3h1Gfjc+8C3pX9n/RoWKQmv7zG4QBEFoS00WXBWEdiQCnm6i/hvCkIt+jUqWuDM3G62kYr81lwOlv1DubHxRUUEQhK7EaDSiKAqKooieaKHDiICnG4pNGklp6kiSfD5+5w4EQh8d+A9qxLckQRAEQWiMCHi6oajYfjiGXIxfVnFj/kEiZS0FznK2F27s7KYJgiB0KYqi4Ff8+Pw+vH4vHp8Ht8+N0+vE4XXgV5qvFSP0HCLg6Y6iM0hIHE5x37MwKQrX1zgA+OzwquAhdo+9s1onCEIP1p3muTi9TuxeO3aPPfivw+vA6XXi9rvx+Dwi4Okm2uJ1JwKe7igihojY/jBoNm6NgatK8kmR9Vi8dcWj8q353eqNSRCErq22kq7d3n2+THkVL4qiIEkSKkmFWlajUWkCNznwfMT7ZPdQ+7o7vqJza4ilJbqr6HQS4gZRNOQ80nZ8yh1l5cw31pURL7IVkeHKIFYf24mNFAShp1CpVERHR4fUcGnXla1PkMdbt+yBz+1DkZsOaLx+L3jBp2r/FdqFE6MoCna7nZKSEqKjo09qRp8IeLorrRF1XH/0/WZgP/QDF1jKeMM0gN3Hdiso5FbnEqmNRC2L/2ZBEE5ecnIyQLOLQnY2j89DWVkZAGqdGllqOJDhV/wU5RcBkNIrBa1a26FtFFovOjo6+Po7UeKTsJPVTs88IZG9iIsZwNHhvyL9p7e4szif/x7bVWMtRqvXUmwvppepV7OnEQRBCIckSaSkpJCYmNjoApKdzeFxkJWfxU033QTAx998jD6iYZFEj8MTPOaD1R8wNm1sl+ytEgI0Gk2b1GoSAU93ptIgxfUjOm0yloPfcXph3UJ27+15nwXTHiG3OpcYXQwGjaETGyoIQk+iUqm6ZLHAAlcBVsXKkSNHAhvUIGsbSVX1EDzGiRO1NpDbI/RsImm5uzMmYo4fgHX4ZSGbc+xF/JT7DXaPnbyaPJGYJwhCj2bz2Ci0FmJU1xUy3Lpxa6MrdNfnVby4fK72bp7QBYiAp7uTZYjuTXzqWMp6jQnZtfTg/5A8Dgptha2qwmyz2ZAkCUmSsNlsLT9AEAShkxVZi1i9YjW/Pfu3wW03X3kz5405j9UrVjf5OI/fg9vv7ogmCp1MBDw9gSEWXWx//EMuCm7qLxuw+118tPdDJAXyqvPw+LvemLsgCMLJqnHXsHTZUhbevJCSwtCE6pLCEu7+v7ubDnoUcPtEwHMqEAFPTxGVTmz8kODdBwpykZFYX76D/LJdVDgrKLIWdWIDBUEQ2sfR6qO89PBL0MjIfe1w/tP3P93o8JYkSbi8YkjrVCACnp5CZ0Id2y94d7jTweXuwH/v29lL0fr85NbkYnVbO6uFgiAIba7aXc1X335FWWFZk8coikJRfhFbN25tsE8tq7F6xfviqaDTA561a9dy8cUXk5qaiiRJfPrpp8F9Ho+H+fPnc9ppp2E0GklNTeXqq6+moKAg5Bx9+vQJ5pzU3p566qkOfiZdQGRq8EevSsfdBYeJlTQUuav4+tD/cHrsHLUeFQnMgiD0CIqiUGAtoLSwNKzjy4obBkUaWYPD4xDvi6eATg94bDYbo0aN4rXXXmuwz2638/PPP/Pggw/y888/s3z5crKzs/nVr37V4NhHH32UwsLC4O3222/viOZ3KcaoGBRLPs6dy6kcdxkmRWHBsQJhnxV8j8taQqG1dQnMgiAIXVW1u5oSewnpaelhHR+fFN9gm0bW4PK5ROLyKaDT6/DMnj2b2bNnN7ovKiqK1atDE81effVVJk6cSG5uLhkZGcHtZrP5pKsw9gimZHSxAzD1OZuK/G2cX7qP/8TE8qPaz3vZH3Lz2Ds5YjlClDZK1J0QBKHbqu3d8fq9TDprEkmpSZQUljTaUyNJEkmpSYybNK7BPo2sweV34fa50al0HdF0oZN0eg9Pa1ksFiRJIjo6OmT7U089RVxcHGPGjOHZZ5/F6/U2ex6Xy0V1dXXIrSurrcisKApGo7HpA2UZ4voRZ+5F9djf41PreKjwKFpkdllz2Z2/gUpnOYW2wo5rvCAIQhurclVRYi8hRh+DSqXivifvC+w4rmBybQXl+U/Mb7RYokpW4fWLWjyngm4V8DidTubPn89vf/tbIiMjg9vvuOMOli5dyrfffsuNN97Ik08+yb333tvsuRYtWkRUVFTwlp4eXpdot6A1IiUMIjmqN0dHXEq618eNVVUAvLf3f5zb93x6R/Umvzy/c9spCILQCvVrhB0oPoBf8aNVBdbBOveic3nh7RdITE4MeUxSahIvvP0C5150btMnlgJrcAk9W6cPaYXL4/Hwm9/8BkVRWLJkSci+u+++O/jzyJEj0Wq13HjjjSxatAidrvEuygULFoQ8rrq6umcFPeYU9LEDMPaeQkX+Nv6vdB//M0dyyO8MHnLQcpCYSLHshCAI3U+po5SU2JSQbededC5nTD2DSf0mAbBk6RImT5/c4jIYEhIOn6Pd2ip0Dd2ih6c22Dly5AirV68O6d1pzOmnn47X6+Xw4cNNHqPT6YiMjAy59SiSBLH9iI/qTdWY3yGpdTxUUhxySLXTwoGqA+KbjSAI3Y9EsHenvvrBzbhJ48Ja80sja7C77W3aPKHr6fIBT22ws3//ftasWUNcXFyLj8nKykKWZRITE1s8tkfTGpASBpES3Y+8EZcyweniAmvdH3Wk10OJvYRDlkP4FX8nNlQQBKF1YrQxbXYujUqDzWsT74M9XKcHPFarlaysLLKysgDIyckhKyuL3NxcPB4Pc+bMYcuWLXzwwQf4fD6KioooKirC7Q5MIdy4cSMvvfQS27dv59ChQ3zwwQfcdddd/P73vycmpu3+ILotUxIRcQMDQ1sJg7ijojK4641dbxODmryaPI7WHO3ERgqCILSOWtVyRkapoxSvv/kJLBDo4fH6vWKJiR6u03N4tmzZwvTp04P3a/Nq5s6dyyOPPMJ///tfAEaPHh3yuG+//ZZp06ah0+lYunQpjzzyCC6Xi759+3LXXXeF5Oec0iQJYvuSYC/n0KjfEl+8KLjr55pD/Hv3+/z6tGvJseSgU+lIMiZ1eBNtNhsmkwkIBMDNzkITBEFoRv0ZqH/+/s/IOpkoXRRx+jhi9bHE6mMxYQp5jEbWUOOvweVzoVfrO7rJQgfp9IBn2rRpzVa4bKn65dixY/nxxx/bulk9iyYCKWEwqY4KDg67EHg3uOuLkp9IPJDApAEXcaDqAHq1HnXnvywEQRAaCKe3ZvXhutptKlmFgoLFZcHisnDIcggAn6NuTa0f1v3AjHNn4PP7xALLPZz4ZDtVGBOIiB+MIeNMagOem2rc/F2n570jnxNvSKRPylj2V+6nt75357ZVEAShEQ5v8zOprG4rP+T/ELz/j3P/gU/jo9xZToWzggpnBetXreeLF74IHnP37+8mKTWJmx+6mSFXDWnstEIP0ek5PEIHqR3aiu4b3PR/5cVc6tWgAK/sfZ/K6jwsLgs5lpzOa6cgCEITnF5ns/u/zv06ZIkISZKI1EXSN6ov45LGwQ74zwP/wVZuC3lcSWEJj9z8CJ8t/6xd2i10DSLgOZWodcjxA4N3/SoND+cdZKJsxq14eT7rVfxeB8X24mZOElC/AJjNZmvxeEEQhJPVXMDj9Xv58vCXTe73+Xw89ZenGk2TqN326IJH8fl8DfYLPYMIeE41xrrF8w6MzEQNvHxoN/01UVi8dp7b8jyqejMz165dK94ABEHoEiwuS5P7NhVuosJZQaS28ZpqWzdupbigmS9zChQXFPPd99+dZCuFrkoEPKcwY9pE8vqehVFRWHLkEPFqI3t+OMjlU38TPOaCCy6gT58+LF++vBNbKgjCqc7j92D3Nl4cUFEUVh5aCcCMjBmNHlNWXBbWdY4WiBIdPZUIeE5hqaZe2E/LpCq2DykuG5d8c5S8V/NwVoR2G+fn5zNnzhwR9AiC0GmcXmeTC3zuq9zHQctBNLKGczLOafSY+KT4RrcfLy6x5eK2QvckAp5TmJwwhL4RSeRPvA67xsTzywoaPa52fHvevHlieEsQhE7h9DrxKo1PS1+ZE+jdOavXWZi15kaPGTdpHEmpScHV048nSRIJKQmMmzSubRosdDki4DmVmZPQJg2ntymVj/TncLS6+XpIeXl5rFu3rgMbKAiCEGD32KGRt6gSewk/Ff4EwAV9L2jy8SqVivuevC9w5/iY59j9Wx6+BR/iS11PJQKeU110b0zxQ3B7osI6vKCg8V4gQRCE9mRxW9CqGy4W+uXhL1FQOC3+NNIj05s9x7kXncsLb79AYnLoOouRCZG88PYLnH3B2U3mCQndnwh4TnWyDPEDGdRvWFiH62NE2XVBEDqWx+/B7rGjU+lCtju8Dr7J/QZovnenvnMvOpfP1tfV2+l9d28ue/syzr3oXLSyNtCTJPRIIuARQK1jyoVXkpac0KCnt77ElETih8dT6axs5ihBEIS25fA6cPlc6OTQgOe7vO9weB2kGlMZlTgq7POpVKrgz8bBRoocRUBg1XSn1xnWEhZC9yMCnlNM/aTj+jV2VMZYFj+3qOHYdj2jrh2GDx8Hqg6Ib0GCIHSYYBBSL4dn84bNrDwYSFae3W82shT4ODMYDews28nOsp0YjIawzl9kK8Lr96KRNXj8HrFqeg8lAp5TyPLlyxk2rG7o6vgaO5lXXceyt/9KSmJsyOMSYvX0vq03hQPL2HjkG6rd1RyqOtShbRcE4dTl8Dj4YdUPXHLmJcFtt/72VtbdvA5Xlospvaac8Lm1shaf4qPEXoJaVuPxiYCnpxIBzyli+fLlzJkzh/z8/JDtx9fYyfzDDexd97/g/pW/i6DwVg33nzUSgKWHPiWvPJsSR0nHNV4QhFPasuXLeOTmRygpDH3f8VZ62f/SftatOvHZoymmFADyrfnIkowiKU3W+xG6NxHwnAJ8Ph933nlns2vIBGvsyCpUiYOD+3vN+hUqWeL6Pd9zYUQ6CvDazjdxOEQejyAI7c/pdvLM/c80OiUdAAmevv/pE64RlmpMBaDAemwGqkLIAqRCzyECnlPAunXrOHq06XLpDWrsaOpmYmkGnseR/mcD8OieHzlNF4/T7+avv7wePKbUXnpS7Wsqr0gQBOHr779u0LMTQoGi/CK2btx6Quev38MDoFapsbtFjmJPJAKeU0BhYeEJH5emj8Mz8tcUpo1Fq/h45eAuEtVGSt11i/gdqDpAhbPihNrWUl6RIAintryjeWEdF+5aWcer7eEJBjyyGpvXdkLnEro2EfCcAlJSUk74ODlxGH10sVjG/x9liYOJ87h4Nb8AnVQ3rdPndbOvYh9VzqpWtSvcvCJBEE5dUQnhFUUNd62s49X28BRYC1AUBa2sxeVz4fF5Tuh8QtclAp5TwJQpU0hLS2t2DZn09HSmTGlkpkNUKuqkEfTXRlM66QYs0ekMtVt4uLLuG9DW3G9xeu3sq9xHtbs6rDa1Kq9IEIRT1qBxg0hISWj2/Su5V/IJr4GVFJGEhITD66DSVRmYqeX3iDyeHkgEPKcAlUrF4sWLARq8adTef+mll0KKcYWIzkCTNIJ+unjyzrwVuzGe6RV1Y+pvfPUxe3LWYXNVk12RjdVtbbFNrc4rEgThlOP2uXH5Xdzz2D2N7q99/5r/xPyQ9y+v30uZo4xqd3WLRQTVKjVJxiQg0MujkTV4/V4xU6sHEgHPKSIzM5Nly5aRmpoasj0tLY1ly5aRmZnZ9IMlCWL6oIsfSj99In/zTmXokrqkvsMvHOHOy57g06XvU+OoYl/lvhYLE55MXpEgCKcGh9eB2+/m/F+dzwtvv4Ax3hiyPyk1iRfefoFzLzo3ZLvFZcGkNiEjU+WsotBaSKm9FKvbis/fsNe4l6kXEMjjkSQJRVFELZ4eSN3ZDRA6TmZmJjNnziQqKjAmvnLlSs4777yme3bqk2WIH8Cq/33Bn+5/h+NHojyVHl5e8B5aWc2sX1/JPmkfg2MHE6GOaPR0J5NXJAjCqcHpdeJTfKhlNededC6f8infX/M9AEuWLmHy9MkN3r9qg5VBMYOIj4jH7rVjc9uwuCxYXBbKHeVYbaG90KmmVLYWbw1OTZckCZdX9PD0NCLgOcXUf3OYOnVqeMHOMT4F7nzkuQbBTn0vP/0+Z08fTzkK+yWZwbGDGyz4B3V5Rfn5+Y3m8UiSRFpaWuN5RYIg9Cg2mw2TyQSA1WrFaAz05Ni8NqRj692UOcooddWVwBg3aVyj7192rx2D2kC0Lhq1rCZSG0mkNpIUUwoevweb20appe48Dq8jpIcHQCWrxKrpPZAY0hLCFsi7yW/2GHeFh798+gpU51FiLWR/5f5GZzucdF6RIAg9XrWrGq1KC8Du8t1hPcbmtpFoSESv1jfYp5E1ROujSTXVDe27fe5gwFPbw6OVtdg8tka/jAndlwh4hLCFm09TWVHDU7vfRltTTFHNUQ5UHcDjbxj0nFRekSAIPZrb58busQd7iMMJeGoTlOMi4sK/jt8dDIAqnBXYPXbUshq3z93o+5bQfYmARwhbuPk0SWY1RW4Lz+55G6OtnPzqXA5UHsDhdTQ4NjMzk927697IVq5cSU5Ojgh2BOEU5/A6cPlcrQp4atw1ROmiiNRGhn0dv9+PUWMkWhcNQKGtEK1Ki8fvETO1ehgR8Ahha6meD0BalIqPTBZi/AqHnWW8uOc9Im2VFFTnsqtsV6MVmU8mr0gQhJ7J4XXgV/yoZBWl9lJK7CXIUvMfWU6vkyRjEiq5+fcQo9GIoiiU2cswGA0oihLs5cm35qOW1XgVr5ip1cOIgEcIW/N5N4HbbfdcQS+dmb8VFmH0K+yxF7Bk34fEOSzYnJXsLNtJXnVeo1NDBUEQatm99uD7TG3vTt/Ivk0e7/A6iFBHEKOLCfsaOpUuENz4vQ0SlwER8PQwIuARWqXJvJteaSx7/Wn+cNH5bDvzZvqpzbxSXIJWUdhSk8Pivf/EUF2Izu8nuzKbfZX7Gh3iEgRBgEDCcu1w1q7yXQAMiRvS5PFWt5W4iDgMGkPY19CpdWhlbUgeT0FNQXC/GNLqWUTAI7Rao3k3hw+T+X+3kxo/nAxjCj+feROnqcw8W1KG5ljQ88CuN6gq3E68rCXfms/Osp2UO8o78ZkIgtAVuXwubB5bgxlaQ2IbD3h8fh9+xU98ROvW09LIGvRqfchMrXxboIfH5/TRN7ovkiRhs4nFRHsCEfAIDdSObyuKEqyHcbxG8240EZA8gtSkkWQYktk6+UYmqSJ5t6CYJJ+PAncV9+97l62HVpGMFrvXzq7yXeRVh7casiAIpwan14nb50av1lNiL6HMUYZKUjEoelCjx1s9VsxaczDxuDXMGjMevyfYw1NsK8br96JWiTJ1PY0IeIS2pdZB4nBSE0fROyKRnyddT19jMv8+WsgEpwuX38NLh//HhzvfIsbjQq/Ss79qf4untdlsSJLU7LetcI4RBKHrq01YliU52LvTP7o/OnXDIqYAdo+dJEMSarn1QYpBY8Cv+InTx6FT6fApPortxWgl7Uk9B6HrEQGP0PbUWqTk4fRKHkuGPp6tk65HSh7B3wuL+b+qwGrq/yvbypNbX8RblUuCrq5mRmOzuARBOLXY3LZgwvKuskD+zrC4YY0e6/Q60aq0xOjDT1auT6fWBas51y9AeCLBk9C1iYBHaB8qDVLSMGKihjPt9DuJu/EHDqdP5e7KKp4rLiUCiV22oyz4+XmOFGwKPmxvxV4KrAWiwqkgnKIURcHitqBT61AUJdjDMzxueKPHWz1W4vRxmDSmE7pe7Uyt+sNa+db8Fqe2C92PCHiE9qNSIyXWJRkeGjyDA6PmcJ7Dyb+O5pPhl6jwWnliz1t1D5FV7K3Yy8Gqg6LKqSCcgtx+Nw6vA51KR7G9mHJneSB/J7Zh/o5f8eP1eUkwJDRbH6w5OpUOnaxrdIkJoWcRAY/Qvup9SxqsiaU0fRy/nHEdfRQV/87NZZrLj0/xB4/Re3xE6aLIseSQXZGN3SMW8BOEU4nDU1dhubZ3Z0DMgEYXIbZ5bJg0plbV3jmeWlaj1+hDpqbn1zS/ZqDQPYmAR+gw8RmTGWZMxx3Tl61n3YZOH83LBUe5qbouqHnipyeprjxCoiGRYlsxu8obr84sCELP5PK5UBQlJGG5qfwdm8dGojERjUpzUtc0a814fJ66Hh6bGFbviUTAc4oJZ8p5u4lMIar3mQyLGYzWEMums27GEZ3BtRV1tXgOO0v4y5ZF7DryHUn6+MDU9bJdFFrDW7i0I4lZYYLQ9mxuG7IsN5q/YzAa2Fm2k51lO1Hr1aglNbH62JO+pkEdWF4i2ZiMLMk4vA6qXFXB/f56vdBC99XpAc/atWu5+OKLSU1NRZIkPv3005D9iqLw0EMPkZKSQkREBDNnzmT//tBpzBUVFVx11VVERkYSHR3Nddddh9Vq7cBnIYTNEIsxYxJDEscQrTHx46RrsSSPCO4eKOmw+pw8tfsffLLzLWJVBjQqDdmV2Z3YaEEQOorFbUGr0lJkK6LCWYFaVjMopmH+To27hmh9dKsWCm2KTqVDlmRUkookQxIQ6OWpJSou9wydHvDYbDZGjRrFa6+91uj+Z555hpdffpnXX3+dTZs2YTQaOf/883E6ncFjrrrqKnbt2sXq1atZsWIFa9eu5YYbbuiopyC0ls6EPm08g3tNJlU28POoy4K73jp8gEvcMgrw0dFveH7Tk8hOa0hBsRp3Tce3WRCEduPz1a2tt2nDJjRogr07A6MHBisu11IUBbfPTZIh6YSTletrbKZW/cRlsaZWz9DpAc/s2bN5/PHHueyyyxrsUxSFl156iQceeIBLLrmEkSNH8t5771FQUBDsCdqzZw+rVq3iH//4B6effjpnnXUWr7zyCkuXLqWgQGTad1lqHZqUkQzsPZU+6rrppGp1BI/nH+ahKhsaZLZaDnD/jwspKd0TPGZX2S4KrAWim1kQeoDly5czbFhdjs49f7iHX034FSs+WwE0nr9j99oxqA0nVFm5MTqVDq1KGzJTq8haFNwvAp6eodMDnubk5ORQVFTEzJkzg9uioqI4/fTT2bhxIwAbN24kOjqa8ePHB4+ZOXMmsiyzadOmBues5XK5qK6uDrkJHUxWoUoYwrBhF5C7/UO++/EFss+7B0d0Or+uLOe9ggISJS1F7ioWbn85+DBJkthTvkfM4mqEyCsSupPly5czZ84c8vNDZ0WVFJaw+tHVWLZYGg14bG4biYZE9Gp9m7RDJaswqA24ffUWEa03pCVKZPQMXTrgKSoKRNhJSUkh25OSkoL7ioqKSExMDNmvVquJjY0NHtOYRYsWERUVFbylp6e3ceuFsEgSckxv0vvNZHjUABS1jh/PvJHyvmcywuVm2eFDTPSp8Sh1Xd46r4e4iDgKbAX8UvoLRbYi0dsjCN2Mz+fjzjvvbHQ2VO22og+L6BfZL2Sf1+8FIC4irsHjToZZG1hTq7aHp/5ECafH2dTDhG6kSwc87WnBggVYLJbgLS9PLGDZqUwJxPedxojk8SSjIWv4BRwadxVRkoq/5x7iamtdl/KDGx/lUOFmkvQJ+PCxu3w3+yr34fA6Oq/9giC0yrp16zh69Gizx3gqPOz4aUfINqvbSpQuqk2SlevTq/UoihLs4alyVwX32byit7Qn6NIBT3JyMgDFxcUh24uLi4P7kpOTKSkpCdnv9XqpqKgIHtMYnU5HZGRkyE3oZDozhrQJDOo9jUHqaAqShrJ1yu14TYncVlb3Gij2VPFI1su8v+01IhSJGH0M+TX5wTV3BEHo+goLwys1UVZcFnLf4XWQaEhs86UfdKrAmloR6ogGuUFOn1P0IvcAXTrg6du3L8nJyXz99dfBbdXV1WzatIlJkyYBMGnSJKqqqti6dWvwmG+++Qa/38/pp5/e4W0+VbRbPR+VBlXCYNL7z2REVD9Uxjg2nHUzVSkjg4ec59OgACuLNjB/3X3kFGWRZEjC5a+bOur0ii5oQejKUlJSwjouPik++LPb50Yja4jUtf0XVJ1Kh1atDQ5rKf66obatG7ficIse5O6u0wMeq9VKVlYWWVlZQCBROSsri9zcXCRJYt68eTz++OP897//ZceOHVx99dWkpqZy6aWXAjB06FBmzZrF9ddfz08//cT69eu57bbbuPLKK0lNTe28JyacOEkCczLRvacwPHEcvdQmto/ODO5+PO8QL5fVECfrKXJV8si2F/jnttcwSXWl53eX76bMUSaqpQpCFzVlyhTS0tKanVaelJrEuEnjgvftHjtmrfmEFwptjk6lQyNrcPvc2H62sf8vdfXe7vnDPQzqP4jly5e3+XWFjtPpAc+WLVsYM2YMY8aMAeDuu+9mzJgxPPTQQwDce++93H777dxwww1MmDABq9XKqlWr0OvrsvM/+OADhgwZwowZM7jgggs466yz+Pvf/94pz0doQ/pIdGnjGZhxNgPlqODmNWUxTLVU8FnOfmYpBhTg88IfeGDtX4LH/Lj+R7KKsthftV/09pwAMdtLaG8qlYrFixcDNBn03PfkfahUdUNXTq+TBEMCstT2H10qWYVRbWTNijWseHgF3kpvyP7CgkLmzJkjgp5uTFLEV2AgMFQWFRWFxWIR+TxtyGazYTIFvo1ZrdYTGv5a/p//cMftt5FfWDfrLiU2gldmSlw+VM3X0YnckaOw7/0jIW9SiSmJ3PTwTVx0yUX0iexDfER8mxQpq9UWz62ttVWbuuJzE3qm5cuXc/vtt4fUTdPEasi8J5MHb3gwuM3j82BxWRidOJooXVRjpzpphyoOccaIMygtLG10vyRJpKWlkZOTExKICZ0r3M/vTu/hEXq2+hVU165dG3I/HMuXL2fOr38dEuwAFFY4mPORnX8fisCyMZ89Lx9s8I2spLCEx25+jC//9yU7y3aKmVyC0AVlZmay7ud1wftD7x3KoOcGMefyOSHH2Tw2zFozZq253dqStSmryWAHAtPl8/LyWLduXZPHCF2XCHiEdnN8BdULLriAPn36hN0l3FydDgAJuO1LJ7etVmiqm1JRFF5+eDEmjYm8mjx2lO6gxF4icnsEoQupv1aV1F9Cp9YxIHpAyDFOr5OEiPYZzqpVUVIR1nHhzjATuhYR8AjtoqkKqvn5+WGPg7dUp0MBykosFFY2n6NTWljG4qXPEq2Nxu13s6tsF/sq94kqzYLQRdRfmRxgUOwgNCpN8L7H50Etq9tldlZ96WnhFaANd4aZ0LWIgEdoc+FUUJ03b16Lw1tt+S3q+30b+dN389hbuoMoXRRHa47yS+kvFFgLgpVbBUHoeE6vs0Fhv+Fxw0Pud8RwFsA5Z59DYkpioPu4EZIkkZ6ezpQpU9q1HUL7EAGP0OZa7JkJcxy8Lb9FxcdFUeG2sHj7X3n2xyfx+jwoKOwu382usl1UOCtaPcwVTn6SmO0kCM2ze+0NFuc8fv2sjhjOAtCoNSxYtKDxnceCoJdeekkkLHdTIuAR2ly4PTMtHddSnY7aGRPNHgOkR0p8py/nqoi+aCQVO6uymf/DfD7LXoZJY6LCVcGO0h3sr9of9jDXyeYnCYIQYPfYQ75saGUt/aP7B+931HBWrczLMnlkySNEJ0aHbE9ITuCtD94iMzOz8QcKXZ4IeIQ2F27PTEvHNVeno/b+4sWLm63loQBPX96LKK+d+3Z/z1KrmvHGdHyKn89zv+Ke7+5iZ8l2DBoDedV5YQ1ztUV+kiAIARaXBbVKHbw/MGYgarnuvs1jw6QxtftwVi2dSseUWVN4c82bwW1//ddf+ef6f3L+xed3SBuE9iECHqHNhdMzE+44eGZmJsuWLWtQNTstLY1ly5aRmZnZ5DG9kuN58embSL71L+SMuAS/Ssugkv28uXsTT+r6kaKNweKx8vedb/L4+oexOKtaHOZqq/wkQRACvTfV7mr06rpCskNih4Qc4/Q6STQktvtwVi29Wo8kSSSb6tZi7D+2PyqVCqdPFDHtzkTAI7S5cHpmWjMOnpmZye7du4P3V65cSU5OTkjXcmPHHNm3mzt++38MVZkp7XcWG6bfTVXycGS/l4v3fscnRwu4LnIYEbKWQ9Y8Hv7xEd7b8SZqWU2lq5IdpTvYV7mPand18LxtlZ/UXk627pEgdCS7147T60RXb1mYgTEDgz939HAWHFtTS6UNWSy0wFqARtZg94qZnd2ZCHiEdhFOz0xr1A+Opk6d2miw1OAYcwJy6liS+0xlpKkPKREJ/DJhLrvG/wFPRDRGWynztq9iWY2KmZEDkYAfijZxz3d380Put2hVWo7WHGV7yXayK7Kpdle3WX5SexB5RUJ3Y/fY8Sk+Ktx19W/6RfUL/mzzduxwFgQCHp1Kh9tfl0hdYCtAo9Jg89hEDa9uTAQ8QrsJp2em3anUEJ2Bvs+Z9O91BiO1sfhSRvHD9HvIGzILv0pDRkk2L2z/hr+RyuCIZFx+D/8+sJwH1i3gUNUB9Go9R2uOklWShcfoCeuyHV2nQ+QVCd1Rtasataxmf2XdQp1alTb4s8vjare1s5oiSRJGjTEk4Cm0FqKRNXh8Hjz+8N4DhK5HBDxCuwqnZ6ZDaI2QNJzovmczPG4YQzRRFA6ayYZz/kx52jgkFCbl/MjS/TtYEDGAWI2JEmc5L2e9ysL1D5JfcxSD2kDKiBQSUhLaJD+prYi8IqE78vq9VLmr0Kv1IQFP/f0qWdVu62Y1x6Q14fXVTVwosAWGtLx+b0hVaKF7EQGPcOqQJDAlokqbQErGWYwyZZCmT2Dv2N+yZcqt2GIyUHsc/G73N3xaUMZV5sFEyFqOWPN5dutzPLHxUcrd5dzz2D0oKA2Kk51IflJb6Op5RYLQGIfXgdPrRK/Ss69yX4P9Vo+1w4ezaulUupD7BdYC1LIar+JtUDNI6D5EwCOcetRaiOuHrvdk+vSayChdAuaY/myZcht7x/wWjz6KGGsJ9/2ymk+rZS6PGopOUnOgOocnNj3B1oRN/PnlP5OQlBBy2l69ep1QftLJ6sp5ReESBRpPPTaPDa/fi9vvpsBa0GB/Zwxn1dKpdKjkui8t1e5qrG4rKIghrW5MBDzCqUtnhqThmPpMYVDCKEZqYvCkj2fDOfeQO/hc/LKa1NJ9PJL1JcucRi6MGopaktlduY8vzF8wffHZwVM9+c6TfLj+Q8bPHN/h3wDbqu5RfWK2l9DerG4rsiyzr3Ifkk5ixkcz2Fm2E4PR0KnDWRCYmq6RNSHb8q35SJKEyyuGtLorEfAIpzZJAmMcUq8xxPSewvCYIQzVxpHbexqqhyuQFlZj9Uj0KdjBoqyvWOqLZ2bUIGQkdlbVJWQPnzgcrUbLvsp9ZJVkUWQr6rA1utqy7hGI2V5C+/MrfqpcVYHhrIrAcNbg2MHB/bXDWSaNqVPap5W16FX6kG351nxUskpMTe/GRMAjCACyCiJTkNMnkpgxmRHGjOCuTWffSVnaWCQUBh/Zwgvbv+V9KY3J5rry9/evu58Pd7+PRqXB4/ewq2wXO8t2UuYoa/emt2XdIzHbS+gIDq8Dh9eBXq0nuzIbgEExg4L7nR4nCYaEkGGljiRJEka1MWRbvjUfrawVU9O7MRHwCKeUFnNF1FqI7Ysm44zgpgxzOofG/4FNZ99FRdIwJMXHyEPreWH3+uAxfvx8c/R77vp2HisP/BeDxoDFbWF32e6G12gHbVH3qD1me4ncnK6ps/9f7B47Lp8LWZI5UHkAqKuw3NnDWbWM2tCApzZx2e1zd1jvrdC2RMAjCI2p92aXmjSS0fpk0mMHcHDyjWw96xYscf2QfXXJi8+qMxgSkYzH7+W/h1dy1zd3sj73O/Saum7x/ZX7qXHXtFuTT7bukZjtJXQUq8eKLMkcqT6C2+/GqDGSakoN7uvM4axa9esBQaCHRyNr8Cpianp3pW75kDp5eXmkp6e3V1sEoWtKGo5WdpNadZSEqiOUJWjYnzAYd85W4FUAzj6wnvMjtKzsM47X1Q5y3ZX8a/8yVuxZGTzNF998weRpk0kzp5FsSsaoMTZxwRMXTt0jm82GyRT4MLFarRiNgXb0hNleQtdz/OvNYDBQ6axEp9KxtXgrEBjOqp2N5fQ4SY9J77ThrFo6dejU9FJ7KX7FH5xZJnQ/rerhGTJkCA899BB2u0jaEk4xETGQPAJN78mkJAxjlCae5OSxwd3V0RmofB4uPvgjn+7fyQIpHra52fznn4PH/GXuX/jNGb/h/Y/eJ6ski0OWQ9g9XedvqT1mewnC8Zw+ZzB/Z2/FXqAuYbmrDGcBxEXGsf7oejYVbCLSFImCQqGtEEVRRC2ebqpVAc/q1av58ssvGThwIO+88047NUkQuihJAkMspIxG0/sMkmLrFjk8NHUeP0++kar4AWj8PvSfb2LX4n14K0PH+ksKS1h400LWrvyeQ1WHyCrJ4ojlCA6vo6OfTQNtPdurrXV23onQNuweO06fE62sDRYcHBwTCHi6ynAWgEbWoFPp8Cpeepl6AYFhLQgsaip0P60KeCZPnsymTZtYtGgRDz74IOPGjRPj+UKHMBqNKIqCoijBIZgTOaZNSBIY4yFlVHDTCE0UvZLHcmjqPH6adBO3f+WjuXkcCxc8yp78LSiKj/1V+8kqySK3OrdTe3zaepV7QWhM7Wu81FFKlasKlaSif3RgxqPT4yQhovNmZ9UnSRImnQm3zx0S8KhV6i7VMyuE74SSlq+++mqys7O58MILmT17NnPmzCEnJ6et2yYIXZtc9+ejTp9IUvxgRmpjOHpUQ0FV813ernIXL/z3bzy07i/sKdiEz+NmX+U+tpVsY0/5HkrtpZ3Sbd7Wq9wLwvEqXBVoZW1wOnrfqL5oVdrgTMDjZ0d1JqPaiN/vp5e5LuDRyBqsHmsnt0w4ESc1S+u8887jj3/8I5988gnDhg3j3nvvxWoVLwShTof1unQ2YzykjEadMRmXK7y1f3QWhRJ3FX/ft5SHNtzPrqPrUTwOiuxF/FL6Cz8X/8zByoNUOivx+Tuu0nGXWOVe6JFcPhc2tw29umHBQbffjValxaA2dGYTQ+jVehSU4Ayy2oDH7XeLJSa6oVbN0nr99dfZvHkzmzdvZs+ePciyzIgRI7jpppsYNWoUS5cuZdiwYSxfvpzx48e3V5sFoWs6luOTMnRiWIcvdpRRYJF5OyqKUreFfxxYxvuHPuP0uNOYmn42htihHKk+Qm5NLpHaSJKMSUTrojFqjE3m2bSVLrPKfTtoapaaEKo9fk+1+TtmnTnYw1Obv+PyuohQR6BX65s7RYfSqXSoZTXJhmQACq2FSEh4/B7cPneD5SeErq1VAc8TTzzB6aefztVXX80ZZ5zBuHHjiIiICO6/4YYbePLJJ7nmmmvYuXNnmzdWELqDKVOnkpaWRn5+fpMVWZPjDJw1NBFjRTlXVlbwcWQkS2PiyPN7WFv6M2tLfyZBG83UlEmcmXEOTpWTvRV70al0mLVmYnQxGLVGjBojEeqIRq8hCF2Nw+tAURQcXgdHawI1nwbFBiosO31O4iLiOmWx0KboVDq0Ki16tR6dSofL56LSWYlKUuH2udultITQflpdh6cl1113HQ8++OAJN0gQuova4brj1Sb/zpkzB0mSQo6p7Zj5y33/x88zRhFVuIPE3d8x9+4dQB4bHoljZWoqX+lkSt1V/OfIF/znyBcMjerP2WnTGJM6EavHSrmjHEmS0Kl0GDVGYvWxgTffbloA9vjFSs8777we1avUU7W2F8jisqBWqdlXuQ8FhWRDMtG6aAD8fj+R2sj2bnKraFVadCodDq+DFGMKh6sPU2ArIM2cJqamd0NtHkonJibyzTfftPVpBaFbaTL5t1cay979O7f/5veM0cRhSh7LrjNvCu4f4fHxWP4Rvs85zJOVVsZLJiRgj+Ugr+96k9u/vo1Xt77E2qPfU2grxK/4sXqs7K/cT1ZpFttKtgXPZXVbu8WaP2Kx0lOH1W1Fr9KTXXFs/axjvTs+vw9Zkrtkb2WkNhK3P3SmFgqi+GA31KoennBIksTZZ5/d1qcVhA7VFj0OmZmZzJw5k6ioQBG1lStX1p3H78Nkr8BUXUBk4f7gYzaecy8Dyn4hLWcjF1dVcHFVBYVqNf9J6c8KnUy+z8auymx2Hct/kJBIN6czKHYQA6MHkqxJDp5re+l2evl6kWBIIFoX3aBUfldQu1jp8YFZ7WKlYmZYz+LwOYhRxwQDnmD+js+FXq3vkgGPUWPEr4TO1FKr1GLV9G6ozQMeQejuli9fzh133BG8f8EFF5CWlsbixYtb/eHbZPKvrAJTAhjjidDGB48ZFJGAbeB5bBowncjCnfQ+8iMpJdnclpfNrcDuqCR+TBnENp2WvY5iij3V5NbkkluTy5oja/C7/MFzVTor0Tq0FNmKMGqMJBoTidXFolK6xlBRS4uVSpLEvHnzuOSSS8TwVg9Ru+jmwaqDQN0MLafPSaQ2sksG5Tq1DgmJVOOxmVo1+ahlUYunOxIBj3BKaannpsN7HCQJIqKDdxMyzqK3uxybrYiK1PHs7zUGX00BvXI2kpK7meGWYoZbilEkmZrUURxKm8E2rZpsRwn7nCUckooY8c4IAB7f9jiDYwczI2MGp8WdRk5VDrlyLjqvronGtL3mft+tWax02rRp7d3UU0pn5UypZBWHLYcbLBjq9rqJMce0+/VPhF6lR6vSkmhIBKDAVoBaUuPyuvD6vahl8THaXXSddHhBaGct5Yq01OMAMG/evJAPizYXnY7U+wxMfaaSkTyG0YZUBkT2oWTgr9A8WIK0sJpCYwaS4icyfxujN73N77cs4082L0+lzODtwdexoG8mE2OHIUsy2RXZ/DXrr9y79l7W5K7B6XVS461bsX1n6U4KrYXYPLY2z/dp6fctFivtHJ2ZM2XUGIPT0WsXDK193Rk0Xaf+Tn16tR6NrCFGH4MsyTi8Dqxua3BqelsRS6e0PxGaCqeEcHpuYmNju0aPQ+1wlykBTWwf4m1lRNTL8zk05Q7yPaUkHN5Iat5WNI5KEvZ8TsKez7EmDiW+zyRGJ02lPGEi31oO8HXFDircFlYcWsGKQysYEjEkeK7v133P6DNHE6GNwKw1k2BICMyUOcnYJ5zft1istON1Rs6U11s3dXDP5j3sMe8BYEhs4HXo8XvQqrRdMn8HQJZkzFozpY5Skg3JFNgKKLYXk2xMxu1zd9lATWhI9PAIPV64PTf5+flhna9Dexx0ZojtC+l1xQxH6hMZFDUA1Zg/sPvCReyeeA3lCYNQkDCV7CH9p7cYvPoxRmR/yx/UCbzS/7f8OeMiRkf2p3pLNZ/e9GnwXPOvns+Vk69k/ZfrqXHXsKd8D9tKtvFLyS/BY1rb8xPu73vy5MlderHSttJVvrl3Rg/m8uXLGTpsaPD+zVfezId/+BDLFguDYurq7+jUui4b8ACYtWa8vrpFRAtsBWKmVjckAh6hxws3V6S0tDSs83VKj4OqrqKrnD6RyLQJpJlSGamNIz1jKq4Z97PnwkUcGTILR0Q0Ko+DmJx19F23mCHfPc/sksOcvstE7qu5DVZwryyuZMEfF7Dw7wvJqc4JJGT66hIysyuyqXJWhR34hPv73rBhQ6sWKz0+76RdhxZ7oNbkTLWF2t6kwoLQLwjuCjd5r+ZxeP1hAJxeJ1G6qC5VcPB4elXDJSYUSRG1eLqZrvsKE4Q2Em6PTEJCQvfocdBHHuv1OQMp/XTM8YNJVUcy1JBK3KirqLr0VQ5Mu4fC3qfjVevQ2suJ3fU5Tz/7YbOnXfvaWhZvWcxtX9/GO7veCW4vshexvXQ7eyv2YnFZWgx8WpObE+5ipaJWz8nryJyp5nqTaj3/4PP4fD58fh+Rmq5VcPB4tXk8ycZA2Yf8mnxUkgqH19HJLRNaQwQ8Qo8Xbo9Mr169WtXj0OlU6kCuT/IIyDgDqddYTKYkUlDTP3E0xsl3UnDZEg6cfh0ratI4Wt18oOKp8KDLUePxe/jF8gsj3hnB2PfG8n3J93j8HgpthWSVZJFdmd1s4NPa3JyWFiut7Sk4fsixNu9EBD3h6cicqZZ6kwCK8ovYsnELKknVpYezoG6mVoIhAYACawEaWYPNI5KLuxMR8Ag93pQpU8LuuQm3x6HL0RogOh16TYD005EShxKpjSRDpad3xtkUp10Y1mmuVI3i2f6/5dfJZ5JuSMbtd/NFzhfc8/09fHbwM+xeO4XWQraXbmdf5T4sLkuDc7Tm912rqXpFXWLmXA9xIv8vJyrcXqKCggJ0Kl2HJ/66vX4sdg+lNS78/paHajUqDRHqCOL0cQBY3BacPmdwarrQPXSLgKdPnz7BpL/6t1tvvRWAadOmNdh30003tXBW4VRRu7YVhNdz01KPQ5cmy2CIhfgBkHEGpE1EEz+QQcmpLT8WGHfkUyZlf8M1fgPPpv+Kv/T5FUNN6fgUH9/kfsP8tfP5d/a/qXRWkl+Tz/bS7WRXZFNqL6XGXYPH52n177s57ZF3cqrmArX2/+Vkfk/RCdFhHRcZH4lBY2jXgoO1wU2RxcnBEis/51bwU045m49UsKvAQo0zvIAlUheJLMnE6mMBKHWU4lHadmq60L66RcCzefNmCgsLg7fVq1cD8Otf/zp4zPXXXx9yzDPPPNNZzRW6oNb23DRZIbk7UWkCQ15Jw5gy50bSUlOa/nYPpEWpODtNISpvMxmb3mTImseZfWgLzxmGsjDjV4wx9UZBYX3Beh5Y/wDv736PYnsxBdYCtpdu5+fin9lSvIWskixGzxjNG/98g+SU5JDrtLanrK3zTsLNBeoqM6vaWkfkTCmKQu9RvYlPiW+2Nym5VzJDxw0lRtf2BQfdXj+HSkODm+15lRwstVLj8KGWZeIMWlxePy5veIGcQW1AUZTgTK0SWwnVNdUYtcYe9zrpqbpFwJOQkEBycnLwtmLFCvr37x+yZpfBYAg5JjKy+SQ4l8tFdXV1yE3o2bpqz01H9DiodAYWv/Iq0Ni3e0CChx/7MwfPe5CjA6bjjIhG5XUSfeRHMjb9g0t/+DuvVNp5JWocZ5j6IAFbSrby6MZHWbLtFX7I+55DlkNYPVasHiuHLYfpfWZvXv/q9eB1/r7076z9ZS1TZk3B4rLg8rlaTIBuy7wTkQsU0JY5U429dqtcVZQ6S7n38XsDOxqPeZj/xHxkldwuw1nF1U6yi2tCgpvUaAMpURHEGrUYdWrUKhlQcHn9LZ4PAnk8kiQFZ2oV2ArwE95jha6hWwQ89bndbv75z39y7bXXhrxxf/DBB8THxzNixAgWLFiA3d78OieLFi0iKioqeEtPT2/vpgtdQEf33BiNRhRFQVEUjEZjg/0dOfuoyW/3qSkse/1p/jjrHAbGDCZm/LVUXvoKB2csIH/AdBwRMcg+N1FHtzIt6xP+tmsDH7hMzNQkIiOxp2o//zn4KS9sfYG7v7ubv6xbwHu732NDwQaKHcXB66SPSeeg5SC/lP3CtpJtbCnaws/FP5Ndns3RmrphK5+/7kO0rfJO2isXqLv2BLVFzlRTr923//U2CgoXXHIBL7z9AonJiSHniU6K5oW3X+Ds2WejUWnaPOBxe/3kVzkwatXHBTcNqSQZuzu8IS2dWodOpSPJkAQEZmpJShPRnNAldbtKy59++ilVVVVcc801wW2/+93v6N27N6mpqfzyyy/Mnz+f7OzsZj80FixYwN133x28X11dLYIeoUN1RtXbJldwl2Vw1SA5qzDWFGN0VELcMOwJw6maALlle4nI3URcwXaMNSWcVrCbF4GjGjVrEnqzwxjJXslLnqeaSlcVW4u3srV4a8hipu/uepe+CX3JMGfQy9SLeEM8Lr8Lq8NKjbVuuYufS34mxZNClC6KCHUEL7z4Alf85gokSQr5XbUmH0is2xWecH9PTzzxBI888kijr90/Xfcnnn/zec7/1fmce9G5nDH1DCb1mwRA77t7s+SWJWREZ1Dtrkan0qFX6dv0OZRZXVTZPSRHtnxerUoOO4dHr9KjU+lIiAjM1Mq35iPJIuDpTrpdwPPmm28ye/bskG+pN9xwQ/Dn0047jZSUFGbMmMHBgwfp379/o+fR6XTodB23iKIg1NeZK4U32culjwzcotLBVQPOKgw1RRicFpJjBlMTM5CKUVeRU3UIw9HNxOfvIM1ylGsKDgbPV2VOYVtSf34xRpItudlXVUzt4Mkvpb+ws3pn8FhZkkkxppBuTidJnVT3/FEosBWQV5OHRtaQMTmDl95+iUULFlFUWBQ8Li0tjZdeeimsoFCs2xWecJ//4sWLm+4FkuDZB59l5oUzUalUIa+3xOGJpEWlAeDyukgxpqCS2+717fUFenf0ahlVGMGIRi3j9Prx+vxN9gLVkiQJs8ZMXERgplaZo6zN158T2le3CniOHDnCmjVrWuzuP/300wE4cOBAkwGPIHSm1vY4dOjq1pIUGvy4rUhOC5HWEiIdlaRGD6Q6MoOyYZkcrCnAmP8z8cW7iS47SHRNIdNrCpkOeHUmimOGknbstH9InkKpyk6uq5I8RzE2r4N8az751vyQnqDHf3yckb1GMix2GANjBoIEI88ZyZKvlnDZaZcB8P5/3ufyiy4nQhte/Raxbld4wn3+FRUVTe9UAjV2tm7cysSzJobsGhAzIFhR2af4MGvNJ9zWxpTb3FRYXSRFhve60KgkXE4fLm/LAQ+AUWvEoDJg0piweqxUOJr5PbRSZ61gfyrpVgHP22+/TWJiIhde2HxNkaysLEC8eQldV2t6HJYvX84dd9wR3HbBBReQlpbG4sWL2z/hWpIC63npzBCVBm4bKqeFGFsZMbZyPJIWiymZqiEXcMRVjb5oJzGFvxBXvBeNy0p07k/BU/3fwa0o6cOwxk/AGRFPpewl12Mjz2vlUFUBu4/1BRXZiig5UsKaI2sASDWmMixuGP0i+gXPFT00mh3lO0gzp5EQkYCm3tIbjanNBcrPz2/0W7kkSaSlpXV+Fe2TZLPZMJlMAFit1kbzxpoTzu8pJiam+YDnmLLiMgD8Sl0wW7t+ll/xIyFhULdd/o7Pr5Bf6UCjUoXVuwOgUcm4fX5cXj/GMDr89Wo9SNDL1IvsymxKHCUn2eqATv0bP4V0m4DH7/fz9ttvM3fuXNTqumYfPHiQDz/8kAsuuIC4uDh++eUX7rrrLqZOncrIkSM7scWC0LRwg/H9+/c3mSvRXnk+zdIaA7fIVPA40TgtxNvLibeV4leZsGfEYO0zlUN+D77SvagObQa+AMBUmo3Rsg8AtyEOa/IwescPxBbTB2tsCm/wDgB3Dvs/DngL2VO5jyPVRyiwFQRmxNTrBcq35hMVGcXu8t1E6aJIM6URHxHfZOBTW4Nmzpw5J5ULdCK60zf3cH5Pd955Jw8//HCL54pPikdRFP65+5/BbWMSxwCB9bN0Kh0RmrarsFxhc1NucxEXTuRyjCxJKBD21HS9So9GFVhiIrsymxL7yQc8nZHLd6rqNrO01qxZQ25uLtdee23Idq1Wy5o1azjvvPMYMmQIf/rTn7j88sv53//+10ktFYSWhTP7KC0tjTfeeKPrVhnW6MGcBEnDoPdk5N6TMaWOI9nUi8G6OAYlj8E86vfBww+PuITyhEH4ZHVgfa9D68j46S0Gr36UjM3vBo8b49cwN3oETw36A2+c8Rj3jLyFC3rPondk7+AxL297mac2PcXBqoM4vU52l+9me+l2imxFePyeRpvbGVW0u8MaYB5f6O+rpd/T/fff3+JrN7lXMuMmjeOjfR/xTd43wX21U7pdPhdGjRGdqm3yKBVFoaDKgYSEJoyhqeO5w52artajk+tmahXawuupbWo2n6gk3rG6TQ/Peeed1+iLIj09ne+//74TWiQIJy6cb9LXX399s9+ku9TMIpUmUOHZEBtY2NRVg9pVTVRpbvCQlEEX4jfOIc/nxFO4HWPBdmKKdxNhr8BUmh08buD3z6GkDcMWNwBVXF/G682MjzkNe9R4vjjWW6RX6cmtyeXVrFdJNiRzcf+LGZUwip1lO4nWRZNmTsOkMRGhjghJim1yllo79Lh05W/u9ZdDyCrJIt2XTrQ+mkhdJBpZ0+zvyeax8cCTD3DT3JsCNXbqPb3a1+78J+bz5ZEv+WT/J41e3+1zE22KbrPnU3lsmYgYQ+srNmtkmRpn40Hy8dSymghNRHCmVqH15JLcxezBjtVtAh5B6Glqv0nfcccdIUXeamcfuVyusM7TWTOLmswXqZ/0rK6roquLH4JRqSHebcOfOhFbr4lYJMivPoo6fytlr+0kuuwQKl81HPmRmCM/AuCITseWMIgSQ13ZiOeHXsda5wG+KFxPkb2IN3a8Qaw+lgv7Xci4pHFYyixoVBp0Kh1R2ijMOjMR6ggi1BHIcl0PQHvVYurMWXgtURSFvJq8ug0SHKk5Qm5NLgaNgURDItG6aFRSXbumTJmC1WuluLqYUnsp/af2Z9Ebi3jxwRcpKawb1klKTWL+E/PRjdbxZtabAFw+8PJgflawDSgYNa3LL2ru+RRUOVBQ0Kpb37ujVctYXb7g/0tLIrWRwZlaRfaiFo5unpg92LFEwCMInai5b9LfffddWOfoNsn5CQMhQg+uGmRXDWZbKWaHBSLVuM1p1Ay5hHzFg6dkN/rCX4gs3oPZkk9EVV7g5q4LHjIOfM/vU/pzYb/f8nXNflaU/UyFs4L3d7/Pp/s/ZXrGdPpG9SXVlIrj2GwwSZLQqXQorrrzWFwW1Hp1mw2t1OrK39xL7CXkVdcFPEaNEYPRgM/vw+a1cajqELIko/bWfTzsLt+NQ3bg8/sw68xE66O5+NKLmXbOtGCNnSVLlzB5+mS2lW7jha0vAHBB3wu4qPdFPMZjwXO5fW40sqbNVkivdngpqXESHXFi63FpVBJunx+3z49O3XLwGaGOIFoXjU6lw6E4TuiatcTswY4lAh5BaERtheSO0FRdnB45s0hWQUR04BadDh4nuGrQOqqIs5US57ZC0ljcSWOwq1QUuWvwFWxDU5CF6uhuIFCgMGnvFxgPSfjUOgbFD+TK2H6sMEr8x5ZDiaea/x78b/CSUdpI+kb1o09UH9LN6cTKscF9WSVZRDkCBQ6jtFEYtUb0Kj2+MJNYm9JVv7lXu6s5WHUQjbphcrdKVhGpjSRSG4nX76XMUhbcV+4sJzkmucEin/Vfu+MmjSO7KpvFPy/Gr/iZmjaV3w/7PU67M+QxLp8rkLDcRgFPocWBx6ugN51YT5lGJWN1e3F5wwt49OpA4nKKMYVD9kPB7fVno4WrR/6Nd2Ei4BGELqozZxZ1GI0+cDMlQFx/cFvBbUVrr0RrLwNFhoyz8PY5m1KHCx4K5LyUp45EazmIxmPHXLQTc9FObgGu1Rr5PDGDjRF6siUvuZ5qLO5qskqzyCrNAgiZ7bUufx0DEgeQZkqjylWFX/EjSzIee11Ox2dffcbsWbOJ1EWGXSSvK35zd/vcHKo6hNPnJEob1eyxallNpLZuPcL4iPgWVzQ/Un2EZ7c/i8fvYXzSeG4ceWOw5k59Lp+LpIikNik4WOP0UFTtJNrQfFmC5mhUMl6vv1WJyxo5MFPrYHFd0c1vvvuGi2Zf1Kq/x1Pib7wLEQGPcErpyJ6bttBSnk9Xn67aqt+3LNfl/kSmgs8L7ppAArStjEh3Xb6E+cy7sRi0OCpz8Bdsw1i8h+jyg+jdNi4/uofLjx1n0xjYntCbHaZo9mhU7PfZyHXV9Vz89+B/kY8GPpSTDcn0j+6PfZudFc+vCB5zVeZVJKQkcPdjd5OZmUmkNjIwDNRMDZmu9s3dr/jJseRQ5igjyZiELMnsLNvZ8gNb4bnNz+GQHQyNHcodY+9oMqDx+DxE6ppf3DlcxdVOnB5/q6aiN0oi7EVEtbKWCHUEFZsr2P/8/uD2Sy++9IRq53T3v/HuRAQ8gtDFtfXMom4T9KnUEBETuEVngKmu2J3WEI9R5YWovnjMGdQMuZg8vwdH+T70xbuIKjtIVMVhjB47kwv2MPnY43xqHUdNvelz7P5kc3+OUEWhu5IiexHZa7PJezXv+JZQWlTKgusX4PK6OGvWWWhUgRwUlavu/+Db779l9vmzg8spdKVv7gXWAo7WHCUuIg5Zkil3lLP6yGqMGiMjE0aSYc4IK2G3OTWeGvon9ufPE/4c0htkMBqCwZVf8VNjq2mT4Sy720tBlZMo/Yn37tSSkXGEuYioJEls+HIDyx9oWFrgRGfgdeTswVOZCHgEoRvo6FXeuyR1vW/xaeNAK4PLisZVQ6y9jFhXDb6EkVhjB1Et+Sn0ulGqjmAo20dk2QGiy3PQeBzEl9RNgX9h1w+oknpTFNOPbXoj1/7rs8avfSxeWXT/In6X8TsGxg6k+Mdi/rmorqjexRdeTFJqEo898xiXZV7GrItn8dHHHzHvznmd+s29wllBjiUHk9aEVqUlqySL17a9Ro0nkA/1wZ4PiNHFcFrCaYxKGMVpCaeFDGfVpygK5c5yDlQe4EDVAfYW7Q3uSzYkc9/E+5pd/dzlc6FTt03+Tkm1E7vbS6/ok6/WrFFJ2Fzh5W35fD6e/MuTje47mRl44m+8/YmARxCEbqHRnilNxLH8n37gceCsKCY6uS8A1u3/Q5sQhSN2MI4hEnn4cFXm4ju6AwgEKrLfi6H8EP3KD5F72Et1ub3ZNrjKXXz7/bessa1ptCeouLCYG/5wA/nWfGZcOIOMMzP49IdPmdB3AgD/W/E/Zs+a3WEfZnaPnYNVB/HjJ0IdwdK9S/n0wKcA9InsQ7Qumt3lu6l0VbL26FrWHl2LhETfqL4MNQ0NnmfFoRXkOnM5UHWAKldVcHv9fKg/T/wzUbrmc4Nc3kDBQb365FZId3p8HK1yYtKdfO8OBPJ47C4ffr+C3MKyFOvWraMwv+lk886unXOyy4v0ZCLgEQShZ6gNfmplnIFGAxpXDZH2CnBWo8SbqNClURvw5Jy/kIiagxhKD5CzPwtoPuABGO/NYOW/fmp857F47NkHnsU+2E5GdAax1M0KMw42csBygGhdNCatCYPacNJDSU3x+r0ctRyl2lWNRqXh8R8fZ0/FHgDO630evx/2e7QqLW6fm+yKbLaXbueX0l/IrcnlkOUQB0oOBM+1bN8yZF0g10mWZDLMGQyIGUC6Lp0/8ScA4vRxLbbJ5XMFqy2fjNIaFzaXl5TIkwucamnVMvZjM7UitM0Ho111Bp7QMhHwCILQM2kiwGgEY3yg+rPXheS2oa8oDh7SN7IP+oR+OPvOIFbeBR8+2OJpJ+/8iWXlzmaPsZfZ+XjVx5iGmkJ6Qdbnr2eIfwgxuhi0Ki1mrZk4fRwmrQmT1oRGPrkei/pLECz/cjmpY1MpdZayZPsSqt3VRKgjuH7k9UxOnRw8TqvSclrCaZyWcBoQGAL7pfQXtuZuDRYMnJA8gaHJQxkQM4C+UX2DdYvstpYDxPoU5eQLDrq9fvIrHRg0qjYLFjUqGY9PwR1GwNMVZ+AJ4REBjyAIpwa1LnBT6uUCpU9ApZEwOi1cdFYsacnx5BeX0VhOtwT0ilKRGGbKyBBvCpoIAznuugDrX9n/Qj4sk2pMZUT8CAbGDKS3uTc6tQ6j2ki0PtDzo1PpgjeNrGnxg11RFD76z0fcdeddwW3X/vpazAlmYq6IIXJ8JL0jezNv7DxSTM1/EMfqY5mWPo2JsRN5m7cBuHX0rRiMJ5cr4/F5gsneJ6PS7sbi9JDcRr07ACpZwqcoxxYRbT7obGkGHhKk9koVtXO6IBHwCIJwSmk0FyiqF6oEP4tfeok5v/0DksRxM6sC/z752HwMEcAnjSet1nd3TQ4TLYMo1A5hKIFZSgMNKeT4SoIrwH915Cs0soYBhgF8POdjAJbvWI7ZbEYja9DIgeUxTBoTeGBIyhAAcsty0eg1uLwurF4rKz5dwX3X3xeyrhVATWkNNa/WcOEjF/LYzY+1WEunvvqzq9pCWxUctLm8SARWOm9bSlhT0+vPwDueJEkoKNz60K1Ue6qJUQWWVqnf87Z27doTnoHVlvk5p2Kujwh4BEEQAGSZzCuuYpkmomFNlNQUXlp4L5nnjMPnspGWFEd+SXmTPUFpkRLnJNlQ5Wehrrckxpu5udiT0vjRaOJHlY+t7grKvVZ2le8KHvPghgdJjE4g1dSLJEMSSYYk4g3xmPym4DG/lP6CRqfB6rVS46zhqfufahDs1G/Q1je2orq5c2f9OH1OEiMSUcsn97FjcXjQNrMiusvr44nP95BbYUchELgqSuDX4z/2g5/Atn7xRh67ZARqlYxKCuTxhCMzM5M3P3iTu+fdTVVJVXB77VpiY2aM4UDVAYbHDWfV/1Zxxx13BI+54IILTqhej3DyRMAjCJ2s29TFOU5bfWvtapqtieL3o/LYWPz8M8z5wx+RJCUk6KntdHj80T9TOTYBpSwbZ/4eYD0AemsRce5i0oFfA35JJjs6mdW6SO48dg4FhVJHGaWOugKJEDoj6okfn8AqWVFQsO6xUlVc1fQTUqAov4itG7cy8ayJJ/W7OV64vUAen6dNCg56fH5sLi86TdOvs3X7y9iWVxXW+XYWVJNTZmNgkhmtSqbGGV7AA3D55ZcTPTyazFGBoOWFf77AjHNnoFKpUBSFYlsxr//zde657p4Gf98nWq9HODki4BEEodWWL1/eo7+1NlkTRZZBZybzqmtZFhHdsCcoJZmXHrqDzBmng9cB8YOp6X0+tQFPyVl3orfmoKs4hLEyF52rhqGVBWS484MBzxd5hRxNSGK/KYaDWh2HVZDrd1Dqqglex+K2IOtkZElGbw8vl6Ws+LgASvHjV/wn3ePSFL/ix+q2YvfYUctq4iPiidZFn9Q57W4fTq+fWEPTbf5yV6Ai92VjenH2oAQCs8wlZCkw5CQRCExfWrOf7OIa8iodDEwyo1HLOL1+PD4/mmZ6kGpFqCKIMtRNw08ckRh8nUiSRKwulkX3L2r0y8zJ1OsRTpwIeARBaJXly5czZ86cU/5ba5M9QbIMXhd47OC2IdebFZaYPBqDbgJOyY8VKHNU4KrYjyd/HxCo3Bvn85JRUcDkioKQ6x3VxpF+7OdnosYRn9wfgzmdLd6DXEvLOUXxSfH4FT82jw27JzC7SiWp8PoDvRoaVSBfSKvShpUo3RhFUXB4HVjdgd4no8ZI/+j+xOhjMGvNja6t1RpOjw+fv+mA5Ei5jb1FNahkictG9yLG2HTOUr8EI9nFNRytDPwuAj08Htze8AIejUpDhKouHynfms9oRgfvZ23KoqywrJFHBnR2vZ5TkQh4BEEIm8/n48477xTfWo9psieodlFUQyxo6urwkDYeSSsR4awhwlkJWhMYk6iIHkltwHNg1qPoHQXoq/IwVOViqjpKhKOKGGvdh+e07Z9i3CPh1ZrIMKfyQIyewkpno2k8kgQJyQmkDU+izF6GQW0gw5xBtD4anUqH0+vE6XVS7anG6rZS467B4/OAFFhEVCtrj/WMSMF/gWDwUrt0htVjxePzEKGOINWUSlxEHFG6qJOeal+f3eUFmg7Eant3JvaJbTbYAUiLCcw6yzsW8GhUEm6fH5fXT7hLc5k0dXlVH2d/TEZ8BqfFB6b3H9+j1hRRr6fjiIBHELqBrpLns27dOo4ePdrkfvGttQXG+EBtIABFAY8DPA50VaXBQwZEDyAioT/OVBdOfFRKEkddNdgKsoFnAbCZEjC4y1G7rUSV7+OVcyXmfBQIBeq/SiQCG+655Vz62y1EK2rMmgg0khb8gCxh1sfCsYU+vX4vLp8Lh9eB0+ukxl2D1WPFrwTyhxS/gg9fIBH42JUUFCQkYnQxJBoSidJFtcnSEY2pdnrRNdH74vL6+Ca7BIBZw5NbPFdG7LGAp8IB1M6w4tjU9PAYtXUzm+xeO4s2LeLqYVdzfp/ziU+KD+scklkiuzIbl9fF5sLNqGU1iqbz/9Z7IhHwCIIQNlFltvWaDFYlCbSGwI16eTgZZyBrZQweOwaXFZxV4HZgJSZ4iOPchRTrtUg1BaiqjjJ9YA7/NGfx548OU1BTl9ycFinx0iw9md4VsHIdRPeGqF4QmRZYkDW2H+jMoDWCPhq1JgK1JgKj2gC6WDAHgovaAMev+AOBjkLofUCv0rdb1WgIJCzXuLzo1I33HK4/UI7N5SPRrGN0RnSL50uPCQRlhRZHSN6OyxN+wFN/mv/klMn8WPEj7+x6h9yaXOaePpek1CRKCkuarNeTkJxA7PBYCq2FqGQVakmNxWXB5/cxJG7ISRdpFEKJgEcQhLB19SqzPWLmmCYCDEaotyQFXhdSvV6g+Oi+gdXiddEQOxCUs/ndeDUX3eoiauxlAHz+8K84v7cPVVUOWPLBaYGiXwK3IAnMyRCVBpG9Av9GpUNM78Bwm84M+kgkdQSSWoesiQC1HlpYb6o9ODw+XF4fMRGND1V9tTswnHXesKSwavTEGrVEaFQ4PD4KLU4yYg1oZBmrK/yZWvXXBLt+5PX0L+7Ph3s+5JvcbyiwFnDnwju5/4b7G3S91QaGf1n0F3pF9go5p0lrotRWSnZFNkNihzS7GGtjuuLfQFep+SMCHkEQwtZSlVlJkkhLS+uUKrM9euaYWgcRdT08pI4GvS4wE8zjDCRIu6yofHU9a2df/AdUBj2otYEP25pCsORBxSGoyIGKg4EgqKYwcGNz3fklVV1PUGTqsWAoHWIyAr1BOjPoIgPBj0Yf+Fetr5uX3w6cbh8er9JoQnFehZ1dBdXIEswcmhTW+SRJIj02gn3FVvIq7GTEGtCqZawuXzAfrSW1S2zUnu/i/heTZkrj5W0vs7diL2WRZSx4bQH/eOwflBSWBI+trddz7kXnNjinLMkkGBMosZWQXZHN4NjBYQc9PfpvoA2IgEcQhLDVrzJbm6xaq/YD4qWXXurwb5Sn5MwxlRpU5kDwUcvUu+7n1NGglQJBjasmELgY4iBlFEgyqLSBQKmmEKryoPIwVOYEgiGPHapyA7f6JFVdABSZciwQyjjWI2QEfSRozaFBkFofmM5/khweX5PpyrXJyhP6xBJnCjPjmEDi8r5ia4PEZbfP3+TQWX2NzTobkzSGx898nGc3P0uRvYgvzV/y8H8e5tbJtwKwZOkSJk+f3OzfiCzJJBoTKbGVsK9yH4NjB7eYF9UefwNdpWemrYiARxCEVsnMzGTZsmUNa9CkpfHSSy91eGAhZo7VU79XIjKlLkHa7wskSHudgX/dNnBZQNYEApKodMiYdCwQ0oC7JjAMZjl6LBA6dvM6A71ElrzjriuDKSkQDJlTwJwa6CGK6ROYqaY1BXqEgoGQLvCvKvwZXFV2D1p1wwDD7fXzzd5A78n5YSQr15ceE5q4rFHJWI+tmh5OwNOUXuZePH7W4yz+eTE7ynbw2o7XgvtGThyJDx8+nw8FpcHrVqfSIUlSSNBT29PTFPE3EB4R8AiC0GrNViPuYD1h5lg4eRcnlZshq0BnCtxCTuo5Fgi5AsNjblugR0hWg8YIMX2hz1mBnh2Vpm4IzJIPVUcCt4rD4LHVGxo7TkRMvUDo2C06IxAU6YyBQEhrPBYIaY8FQrqQXiGvz4/V5UXXSMCz8VA5NS4v8SYtYzNiGuxvTnpsoNfkaLCHR8br9ePy+EPyyE+ESWvivon38f7u91mZvTK4/frV1yPrmu7xyjBncOOoG+kf3T8Y9BTbipEqJVJUdblx9V8DJ/I30BVzfdqbCHgEQTghTdag6WDtMXOsI8sAhJN30W65GSpN470sXnegN6f25rKCqzrQI6SLDARCnBXoUZI14LaCtQSq8wO9P1VHAsNhtjJwVAZuxbtCryFrAr1Q5uRAEGRKCSRNx/QOBElaY7BXyOFV47I7iWpkxfba4azzhiWjamUydW0Pz9EqB35FCSQ7S+D2tbyIKAReJzWuGn4u+RmNtuHvUSWruGbENSRpkvgzfw7rnLk1uTzwwwNc3P9i5gyag1alJdGQyPLly1nyyJLgcfVfAy6XK6xz1/4NnKq5PiLgEQShW+vqM8eaE07eBdDx+UlqbeDGcWtf+TzHhsWcdcNjTktgiEpnhug0UCYCUiCQ8nnBVgI1Bcd6hfLAkhv42e9pPE8IQB8FpmQwJ4EpGb8uGa23F7qEDBStETQm0OrJr4Ed+ZZAsvLguFY/zaRIPWpZwu31U1LjIjlSj4yMI8xFRAF0ah1aWYvL50LTxBDd2WlnB39eMnMJEYaIkCKOEMiBc3gcvL/7fdYXrOe/B//LluIt3DzqZg6vP8yjNz/a5GvgkUceCautKSkpp2a+2zGS0hWqmXUB1dXVREVFYbFYiIw8uQXuBOFU0FUSGn0+H3369Glx5lhOTk6X6rKvbXdTQxGSJNGrV2DKcnPH1H9unfZ/4vMe6w1y1fUKua2BZGmvG3yuQB4RBIaqJNWx4bGiY8HQ0WPBUB7Yy5u9lFsfhzsiGbchkect0/mgbACnxzl5cHIEis4MagNoIwJJ2bU9WCpdYJiuEbd9+DNHKuw8fNEwxveJpdzqItqgDauWT63tpdupdlcTq49tdL/dZmdi78DCrT8d+QlDIz1V9W0p2sI/dvyDKlcV+OHI/CPUlNY0emz910lLfwMHDhygf//+Hf56au/XZbif36KHRxCEbq2rzhxrSTh5F83trz2mS+QnqdSgaiRHSFHqBUHH/nXbAsGQWh8YuorrHzhOIhCU+H1gKz02RFYI1Uexl+WhqclD47WidZajdZbjqsjmC9cVAPyx+lVGrP4Zjy4WlyEZd0Ri8F+3sRduUy/8uqhAXpJWD/KxHixZTVqMniMVdvIq7YzvE4tGJeNw+/D7FeQwh8hidDGUO8opsZcQqY0Mqc9zIsYnj2dI7BDe2/0eK1evbDLYCfyKA6+T+Q/M55knnmn2b2DDhg3dPt/tZIiARxCEbq+rzRwLR1tWo+6yla0lqW5dseP5vIGen/q9Qi5bIFdIYwjk9SQOw+eHrGJQZBUxshOdswStvZi1hRIVRyJJlKuZqjsAXtC4KtC4KqByd4PLeTWRuA2JgWAoIgl3RBJuYwq96cN6NOQVFEJfP1pFjd0n47IqROgNgeE6uflgOdWUikFjoNhWTLmznCpnFQatAZPGdMILppq0Jm4ZfQv+n/28zustHq9KUPH43x7npYdeorSorkhlr169ePGlF8nMzORf//pXWNfusq+nkyQCHkEQeoSuNHMsHG2ZU9QV85NapFIHbtpGhjdqh8G8Lpx2Oy57JWbJhU9xYNdGYjf15uMDgTyb6f3N7Bv6BiqvA62zBJ29GK2jGK2tAJ2tEK2tALXbgtpTjdpSjcFyIORSZ/gm8SG3U5aXTe/vHselT6BCnYi/KAni08DcCwwxx3qHjgVAKm3gduxntawmPiKeOH0cVo+Vckc5xbZiim3FaFQaNL4TX0B14qCJYQU8/TL6MXrSaMZNGcesobMAePrdpzlj2hnoNDq2l2zHZwpv2Yxu+XoKgwh4BEHoMbrKzLFwhFO1OtzcjNrK1l1lkdmTVps0rTNj97twR2jQRkUEhr18Hgora8gq24MEnDsiDXQefB4HDm0kDlOfQEJ0LUlGVtyBoTBHKVp7CVpH4bFgqJD+tkBvxiFvPObSbURKkABwuF579FGBhV8NCWBKAGNCYFp9ZK/ALLOIaNCakNQ6zGodZk0kqdGRVHrtFDnLKLQVBU9l99qRvTKyFLipJFWzVZ3HTRrX/JpcgCZWw1+r/0rvH3uTqkkNbh84biCKpOBTfNR4akg+LZmElIRAD1Cjy7t1XqX0jiACHkEQhE4QTu7R4sWLAbpdflJbcnp8gFQXFKi0fLUvkNMyJiOapF596w72ecDvDvzr84DPDR4nfo8dpyYSpzEtMJQW/LSXcPlBWqlgwcSOYXeT7MnHX11ArKcYrb0wkHjttARu5QcbaaEUKK5oiA8ERcZAQKQ1p5AUmUKiKZk+OiMHDq6h3O/Ar9hw2pz4JQm/pMInSYEFWKXas0no1XqMGiNqWY1KpeK+J+/j7v+7u8GaXLXSrkrDrbjZX7WfbFd2cPsD6x9A1slISBg1RsxaM33n9qX0qdIG56j9/b744ott/nrqKjV/RMAjCILQScLNPepu+UltqdrhQa2q6wHx+vys2VsMNFJZ+disLL9awer0IqskZD3IkoRK8iP7vaiUY8GQ3wNeFzqPgyRjOUU2hZ0RE5F6eSmxS/QywZA4OTC8Zq8Ae1mgrpC1+NjsssLAvz5XYGaZvRzKsjmehESUIZYoQxz9DPEohjh8xjh8xgR8xiQ8kQl4dVH41Tp8Kg1OFMpc1VQ6qvBJEnqtmWmzp/HC2y+waMGikDW5knslM/+J+cy4cAaFtkKOVB9hf9F+dhPIYdKr9Lhxo6Bg9VixeqwwBNJvS6fwg0K8lXXT72OSYpi3cB7Dpw+nxl2DWWtu8FxORFeq+SOmpR8jpqULQvfXVabKt1bt+w80nXsUzjE9jd+vsOlQBT5FISoikAez4WAZi77YS7RBw9tzJ6BuZDHRSnvgQ14jy/gUBb9fOfYv+PyB3hQJ8CsKiWY9T67cw5Yjldx8VhoXDImmssaBQeVjfJIUmFHmcYDfG+gx8vuOzSqTAktqeGzHAqLyQFBkLakXFBUFHtMSnTnQQ2SIA0M8fkMsDkMMNfooirV6qrQ6vCodHo+K2RPnArDk/ReYPP1MVBpdIBfqmOOnwGsjtFjdVmo8NYF/3TXUuGsorSzl8ZmPA9D77t6YRpgYGj+UaWnTGJkwkl6mXhh8BlLiA/k8J/Kaa6rmT21vUlvV/BHT0gVBELqJcHKPulN+Ultxen24vD6MurqPqtrKyjOHJDUa7ADY3V6GpkSSFmPA51cCN0XB5wv86/X78fshp8yG3e0jPdbAliOV5FX7QB+NRjLhVBQ88YFp6vh9gcDF5z6WUH0sqdrtCKw7ZkwMDYjg2BCVHFiItbYHyF5e10tkLYaa4kDA5KoJ3CpzOPYojMduyYAiq/FGRFMu1QXwp/s3os46hGJMDOQRGWIDU/3d9apEu22odSqiNWai9dGhv6N4O48TCHhmnzObjeUb2Vuxl70Ve+kd2Zv4/fF8/uLnweMvuOACeqX14oUXXyAzMxNZCgyV1c8/8vg9eHwe3D43Do+D226/rUut7yUCHkEQBKFLsrt9uHx+Yo+toVVc7WRbbhUA5w1PavQxTo8PnVomxqhFJUvNLjdRZnVSUu0mPaZ2Ta3AIqJalUy1043b6w8EPLIK5AjQNLFieW1A5HXV5Q753OC2BwIaU+KxafjuYwnVtQk7ciBwclSBoxxstUFRaV1Pkb0cye9FYyvD7K7LvRm0dxVGbb2hPpUWpz6SKqmuFlL05nfQxCXhjYjDY0rEb4gOFGFUaf+/vTsPj7O67wX+ffeZkWa0WNZmywabYFaHshlDWJL4wg2EJCUULiSE0jQ8uTFpg9M2EJo4lBLzpL0tvdSENg3Q55bECbkmTQl1w74kcPPU4LLaibENXiR50TKjWd713D/eZUbLyJItaRZ/P4/9SDOSRkd6X837nXN+5xygUCzs/oOTP4vrlGvw+Hv/jqffexqvP/M6dv/9mA1i4RfPX3vNtVj3j+vw4Ss+HK0WrUgKhBB+4PH8wPPrl36N3n3lp7dXYs0fBh4iIqpKecsNRo/8C/srOw5BAFi+oAldTROHj0zBwbxGA0nj8Je3ZEzD3sFCya7p4SaiEmxXwHQ8NBhTaOhhA5FXDEGlgcjO+z1ADfMB97igl8gGROleXsLv/ckPAof6Adzl391zDoQ9AIwcgGRloLoWGrMHIZWEou43No4ORaoBK5aCZaQwJBeDUePrP0ZLayf+oOFEfOLUU3D5mi9P/HMEnTVrb1uLJ5NPIq7HoTkafnrNTwEAlz5yKfJyHmkrjUMvT75idmgu1/xh4CEioqqUMW1oJT00b+wdBgCcuXjiXdH9XgYX7Slj0qneoZimAJLAwlY/8BzKWshZDhK6CgHAdKa2bs1hyTIgl1mA0W/46CAUFlW7lj9sZmf9hRkbjit+zbn/E1IiFhQjucEmrcPAYB+AewAAVvcZ0J1BKLlDkK0sVMeEOnIAiZED0KziUNOiXz8UBaNn3hPIHByZ9MexDlnY8doONJ7cCM8shrM9I3uineDVpqnFi7lc86fqA8+3vvUt3HnnnaPuW7ZsGbZu3QoAKBQK+OpXv4oNGzbANE1cdtlluP/++9HRMXF3JxERVT/PE0jnHRiaX9/hCYG39qUBAKcvaJrwa7KWi4SuoiWhT+l7xFQFiiQjriloSWgYzNnYM5jHiR1JSABMe4YCz+FIUsmGrWV4HpAeLN5ecCYQU/2aIjsPNGb9IbVUcRsKfeUfQQtDkWv7U+vzQ/Dyg3AG9gH4JwCA3bwYjjsI1cygPz21jVP/YHcaKxbEcEDo+Gxw3zdiJ6CtuR2NifmIfaAL//3hO7F//yFMNDWqEmv+VH3gAYBTTz0VTz31VHRbVYvNvvXWW/Hzn/8cjz76KJqamnDLLbfgqquuwi9/+ctKNJWIKqhuFt4jmI6Hgu2iQfef73cdzGLEdBDXFCyd3zjh14yYNnpaEn7PzRQYmgxdlWE6LnpaEhjMDWP3QA4ndiShyjJGzKnvmj7rZHl0D1GyAxg7C9Fz0eDaEPnh0b1Eru3XEzl5wClAdl2o8zMIA4/20bVQEzEIz8G8rpeBjd85bHMu1Aq45NA+ZEt6iq7c+uyoIbT1Fzv4vR+PXz6oUmtI1UTgUVUVnZ2d4+4fHh7G97//ffzgBz/ARz7yEQDAQw89hJNPPhmvvPIKzjvvvLluKhERzYC87cJyPLQk/CGScDjrlO7UhIXIrudfUtuSUym68RmqDEOVYTkeFrYm8PreYewOC5dVGSOmG80oqgmy4v8vN3QG+D1Fnj26t6hrORDTINl5fPRjC7Cw8yHs7T9QpmcG6GpvwdLf/zPssYYxMrgfwIMAgIMdp8AVI9ALwzAKGXz6ZBU/uSaOP/r3AvZmig9WqTWkaiLw/Pa3v0V3dzdisRhWrlyJdevWYdGiRdi8eTNs28aqVauizz3ppJOwaNEivPzyy5MGHtM0YZpmdDudTs/qz0BERFOXsxx4QkAOwkYYeMoNZ2UKNpIxDc3xqe9bJUkSUnEN+wbzJTO1ioXLpuPBcj0Yah0tASDLgGz4a/+EmhZEvUUKgL9b/0CwujcmXN37vru/gZ6lFwCOiWwmgzDwGOfdAjeuIS1c2MKGU0jjg5cM4t9+rx9n/t79AICNP3oQn/j057jS8kRWrFiBhx9+GMuWLUNvby/uvPNOXHjhhXjzzTfR19cHXdfR3Nw86ms6OjrQ19c38QMG1q1bN642iIioEqYyFHesDdeNFByost+7M9X6nZNa42XX5imn0VDheGKCmVoyRiwHplM9gWcmz4HJHmuqK4BDCCA9FH08edyH0BDXi1t72HnAySM7fwiAH3guXfWxiq0hVfWB52Mf+1j0/vLly7FixQosXrwYP/7xjxGPl5kCOAW333471qxZE91Op9Po6ek5qrYSEdHRE0JgOG/DCNbfOVz9jum40FUJLVOaQz6aockQEFgY9PD0pQuwXX/9HcfxYNoeMMkIUb266qqrsGrVqslX9w6LrUON88fXFQHASMmsL2Nmtqw4EtOLwlWgubkZJ554IrZv347Ozk5YloWhoaFRn9Pf3z9hzU8pwzCQSqVG/SciosoLC5bDnpXD1e9kCg5aEjpSsem/ho9pCjRVRjKmIqEr8ASwb8iv4xEALNeb/AHq2Iyt7l0lNVA1F3hGRkbw7rvvoqurC2eddRY0TcPTTz8dfXzbtm14//33sXLlygq2koiIjlTecv2hJG10wfJEw1lCCJiOi45U7IiKi2OqAkORYbslw1pB4bIiychbVTRTi45K1QeeP/mTP8Hzzz+PXbt24Ve/+hV+93d/F4qi4LrrrkNTUxM+//nPY82aNXj22WexefNm3HTTTVi5ciVnaBER1aic7cKDX7B8uPqdXLD2TmvD1NbeGUtXZcQ0xZ+pFQxrla64PFKYo7V4apTrFn8/L7zwwqjb1abqA8+ePXtw3XXXYdmyZbjmmmswb948vPLKK5g/fz4A4G//9m/x8Y9/HJ/+9Kdx0UUXobOzExs3bqxwq4mI6EhlCw4UaWr1OyOmg/lJfcpr70ykKa7BdD30tIY9PH7g0VUZecuF5x07xeLTsXHjRpxyyinR7csvvxzHHXdc1V6Dq75oecOGDZN+PBaLYf369Vi/fv0ctYiIiGaLEALDhWLB8mT1O64n4AmB+cmjqypOGCo84UVT00tnauWCmVpxvTpmalWLjRs34uqrrx4302vv3r24+uqr8ZOf/GTO19k5nKrv4SEiomOH6XjIWc64wLN8guGsTMFGKj69tXcm4n8vCQua/cCzdygP1xPQgtqeGdtTq064ros//uM/nnBae3jfV77ylaob3mLgISKiqlGwXVi2gKEqcD2BN/f5gee0CQJP1nLRmTKmvfbOWDFNga5IaG0wop3S92cKUGQJrhCwnGN3ptZEXnzxRezZs6fsx4UQ2L17N1588UUAxTV/hBBomGja+hxh4CEioqqRs1x4QkCRJew6lEXWdCes3zmatXfGiqky9CBghb08uwfywUcFTAaeUXp7e2f08+YKAw8REVWNrOmM207i1Anqd45m7Z2xVEVGg67CdFwsDKamh1tMKJJfx3MsKtcz09XVNaWvn+rnzRUGHiIiqhqlKyy/WWb9HSEErKNYe2ciyZgK2yspXB4sFi5nj9HAU86FF16IhQsXlv3dS5KEnp4eXHjhhXPcsskx8BARUVUwHRc5y4WhyZPW7+QsF3FdRUviyNbemUjC8FdZjqamB0NamiKhYHnRbuzkr8D8d3/3dwAwLvSEt++9996K7ZlVDgMPERFVhYLlwQo265ysfmfEdNCW1Gd0qrihKpAALCjZNV2IcKaWx8LlMcINRru7u0fdv3Dhwqqckg4w8BARUZXI2y7coGC5XP1OtPZO48zu6BnTZOiqjPmNBmTJnwE2mLOjqekMPONdddVVePvtt6PbTzzxBHbu3FmVYQdg4CEioioxYtoIs025+p0R00EqpqE5cXRr74zl76mlQAigI+WHqd2DuWhqullla8pUixnbYHQOMPAQEVHFCSEwmLMPu/5O1nQwP2lAO8q1d8aSZQkNMQWW60WbiO4JVlyGAEybPTy1joGHiIgqrmB7yJkOYpPU77iegCxLM967E0rFVD/wtIYztfzCZVmWULDZw1PrGHiIiKjispaDguPB0OSy9TtZ00GDriB1lFtJlBPXVUgQ0Vo80dR0WULWZOCpdQw8RERUcXnLBSAgS1LZ+p2sNTvDWSFDlSFLcrTa8p5gajp3Ta8PDDxERFRxA1kLulK+fscTAgJA8wyuvTNWTFOgqRLak/52FQM5C1nTgabIsFwXlss6nlrGwENERBVlOR4ypoOYVr5+J2e6aDBUpOJHv5VEOYYqw1AVqIqM1gY/WO0ezJXsms7AU8sYeIiIqKLylouC5SKmlq/fGbEctDXoMNTZm/YsSRJSMX9PrXCLiT0D+eLUdId1PLWMgYeIiCoqazlwPQFVkSes3xHCX2ww7HWZTcmYBscrTk0PC5cBcPHBCZTbYLQaMfAQEVFFpfO234tSpn4nZ/lDXLM1O6uUockAJCxsHR14ZHBqeq1j4CEioorxPIGhvD1p/U7WdNDaoCOmzf4qvoaqQJEldDf5qy3vGSxuIsqp6bWNgYeIiComZ7vIWw5iZdbfEULAER7aGo05aU9Mk2GoMjpS/vfrTxdgOR40RUbO5NT0WsbAQ0REFZMzHViugKEqE9bvFGwPMVWZ1dlZpQxVQUzz/zcYCjwB7B3K+1PTPU5Nr2UMPEREVDFZ0wGAsvU7I6aD5oSOhD43gQfwt5iwXVHcU2swB02RODW9xjHwEBFRxRzKWpPun2V7HtqSczOcFWowVLiiGHh2D+SgKjJc1+NMrRrGwENERBVRsF3kLLds/U7BdmEoMprmYHZWqZimQAKwsGX0JqIC4JBWDWPgISKiishZLkzHRUxTsLUvAwA4pTsVfTxrOkglVDTosz87q1RMU6ApErqimVr+1HRFklGwOFOrVjHwEBFRRWRNB54AZEnC9v1+4DmxPRl9vOC6mN8YgyRJ5R5iVhiqDE1V0J7yA8/eoTxcT0BTJIwENUdUexh4iIioIobyFjRZRjpvoz9tAkBUv2M5HnRFnpPFBsfSFBkJXUZTTIuKlfdnClAVGTnLhRCcml6LGHiIiGjOOa6HdN5ff2f7gREAQFdTDI0xfzZW1nSQjGlIGnM3O6tUKqbB8QQWNAd1PAP5IPx4nKlVoxh4iIhozmUtF4Wgfmf7fj/wfKC9ZHd020V70oAsz+1wViihq/Ag0NNanJquKzJsz2Phco1i4CEiojmXsxy4roCmyFHgOSEIPLbrQZWlOZ+dVSqmKZAlYGHYwzPoT013HE5Nr1UMPERENOdGCg4k+L034ZDWCUHBctZ00BhTkYxVMvDI0BQZncFMrd0D/tR0SNw1vVYx8BAR0ZwSQmAgayGmyRjKWTiQMSEBWDq/AYA/Xb0jaUTr8VSCoSrQFRkdyWBq+lAOQghI3DW9ZjHwEBHRnMrbLvJ2UL8T9O50N8eR0FW4noAsSWhK6BVtoyJLaDBUtDbokCUga7oYytnQFJlT02sUAw8REc2prOnCtF0YqjyuYHnEdNBgKEjFKjM7q1QqpgIS0BGsx7N7MOfvms6p6TWJgYeIiOZU1rQBSJAkaVzBcs5yMD9pQFUqf3mKBxuWlm4xEU5N50yt2lP5M4qIiI4pQzkHhupffn5bEnhcz+81aWmo7HBWKKbJkCUpWotnz4Dfw2NxE9GaxMBDRERzxnRcZEwbMU3BoRETA1kLsgQsaWtE3nKRMFSkKjg7q1RMU6CpErqagsAzlIcqS3BcwcUHaxADDxERzZmc6cK0PcQ0Be8GBcsLWxKI6wpGLAdtDTp0tTouTYYqw1AVtCcNAMDugVy0rxd7eGpPdZxVk1i3bh3OOeccJJNJtLe341Of+hS2bds26nMuueQSSJI06v8Xv/jFCrWYiIjKydkuHOFBkaVRw1mAP129ucKzs0pJkoSkoaIt6bfpUNZCznIgATA5Nb3mVH3gef7557F69Wq88sorePLJJ2HbNi699FJks9lRn/eFL3wBvb290f/vfOc7FWoxERGVM5y3oEn+pad0hpbjepBlCXFdqWTzxknFNeiKgpaEP8y2ZzAPVZaRtTg1vdZUft7fYWzatGnU7Ycffhjt7e3YvHkzLrroouj+RCKBzs7OuW4eERFNkesJDGcdxDQFQoiSFZYbUXD8Ya64Vl2BJ6YpAAQWtsQxmLOxeyCHjlQMI6Y/NT0c4qLqV/U9PGMNDw8DAFpbW0fd/8gjj6CtrQ2nnXYabr/9duRyuUkfxzRNpNPpUf+JiGj25CwHeccPPIeyFoZyNmQJOL6tAabtImmoVVO/E4rrCnRVRndYuFwyNd12uRZPLan6Hp5SnufhK1/5Ci644AKcdtpp0f3XX389Fi9ejO7ubrz++uv42te+hm3btmHjxo1lH2vdunW4884756LZREQEf8sI2/Wgq3JUv7N4XgMMVcEhx0JTvPouSbGgcDnaUytYfDBrOTAdt+oCGpVXfWfXJFavXo0333wTL7300qj7b7755uj9008/HV1dXfjoRz+Kd999F0uXLp3wsW6//XasWbMmup1Op9HT0zM7DSciImQKNuRgCChacHB+o79qsQQkjOq7JKmKjAZDxfxgppZfw+NPTedMrdpSfWdXGbfccgsef/xxvPDCC1i4cOGkn7tixQoAwPbt28sGHsMwYBjGjLeTiIjGE0JgMGfDUPwane37MwD8+h3L9aArEhJ6dV6SmhMa2hr960XvcB6OJyAArrZcY6q+L04IgVtuuQWPPfYYnnnmGRx//PGH/ZotW7YAALq6uma5dURENBUF20POKhYsl05JL9geYqpadQXLobimoDmhIa4p8ASwbygPgFPTa03VB57Vq1fjX/7lX/CDH/wAyWQSfX196OvrQz7vn3Dvvvsu7rrrLmzevBm7du3Cz372M3zuc5/DRRddhOXLl1e49UREBMCvebE9GJqM/RkTmYIDVZaiguWmhApFrs4ZTzFdgabKWNBSUrgsy8haDDy1pOoDz3e/+10MDw/jkksuQVdXV/T/Rz/6EQBA13U89dRTuPTSS3HSSSfhq1/9Kj796U/j3/7t3yrcciIiCuUtF0IAcsmGoYvnJaApMmzPQypeHdtJTCSuKTBUBV1B4fKewRw0RULW5Fo8taQ6B0xLCDH5tL+enh48//zzc9QaIiKaLiEEBrJWNKOpuEN6Ep4QkKXqrd8BAE2R0aAr6EiGM7Xy/iaijr+JKGdq1QYeJSIimlWDORsHRkw0BrOwwgUHP9DeGA1zJapsheWxmuIa2lPBnlrB1HTb9Vi4XEMYeIiIaNZ4nsCewRwgUFKwXJyhVXBcJHQVRpX3kiR0FZ2pcEgrD0UGbE5NrynVfYYREVFNG8hZ2J820RJsCtqXLiBrutAUCYtaEzBtD81xreq3aIhrCjpSMaiyBMvxcGjE8qemM/DUDAYeIiKaFZ4nsGfA3+ZnbP3O8W0N0BQZAgKNseqt3wnFdBkJo3TF5WBqusOZWrWCgYeIiGbFwayJAyMmWhv06L7flhQshzukJ7TqDzyG6m9s2pUqbjGhypypVUsYeIiIaMa5nsDugRwUSYamFC81xS0lGoo7pFd5wXLIL1wO6ngGwj212MNTKxh4iIhoxh0cMXFoxEJLori+jifEqCnp1bpDejkNhoqOaKaWPzXddDzYnKlVE2rjLCMioprhuB7eH8hBV2SoJb07+4byyNsudEX2C5Ydryp3SC8nrinobgpXW/Z/PidYi4eqHwMPERHNqP0ZEwNZC80JfdT9Ye/OkvkNkCVU7Q7p5cQ0BYvm+YEnXXCQsxxYLgNPrWDgISKiGWM5fu9OTFXG7Y21vWTD0GrfIX0iMU1GKq6jrdEPcnvDTUQZeGoCAw8REc2Y/ZkChvMWmibYG2vsCsvVvEP6RCRJQlNcixYg3D2QhxBci6dWMPAQEdGMMB0XewbyiGvjdz53PYF3DxQLlgtOde+QXk4ypqJj7NR0i1PTawEDDxERzYj+4QKGyvTu7B3Ko2B7iGkyFjTH4XiiqndILyeule6a7s/U4lo8tYGBh4iIjlrBdrF7MI9GQ4M8wTYR24P9s5bOb4QkAbKEmqrfCcU0BT2tCQD+TC1NlVFwPDicml71GHiIiOio9Q3nkSnYSJXZJiJaYXm+X7+jq3JN1e+E4pqC4+Y1APBno7me509NZ+Cpegw8RER0VHKWgz2DeSSN8puAls7QCndIj2m1dwmSZQkLWmJoDKbT78+YsF0B02bgqXa1d7YREVFV6R0qYMR0kCzTu+N6AjsOZAH4gadWdkgvpymuo7PJX3F572AeAuzhqQUMPEREdEQKtovt+zPYeSiLprheNsC8P5CD5Xr+SsXNcXg1skN6OTFNiaam7xnMQwDIFGwIISrbMJpU7Z5xRERUEZ4nsD9jYtehLIZyFprjOhomWTE5LFg+ob0RniegSLWxQ3o5MU3GwhZ/xeXdgzk0xXTsPJiDocpY2JKo2Z6rele7ZxwREc25dMHGroNZ9A0XoCsyupriE87KKrW9ZDjL3yFdrpkd0icS1xQsavULl3cP5tFgqBAC2NY3AkmSsLAlUeEW0kQYeIiIjkEjpoOhnIW4pqDBUGGo8qQ9E5bjYe9gDu8P5FBwPMxvNKAph6+KeKc3jWe37gcALOvwd0hvTug1s0P6RFRFxgfaGwEAvUN5uJ4/RCcgsLUvA0mSsKA5XuFW0lgMPERExxjH9bC9P4N9w/7CeboqI6GraEloaDQ0xHUFCV2BpsgQQuDAiIldB7M4NOIPX7U2GFP6Pr/pz2Dtz95C3nZxRk8zzj2+FQcyJpoTtX/pWdreAF2RYbkeeofzWNiSQDKmwRPA1r40ZAnoamLoqSa1f9YREdG07BvKoz9TQFcqDlmWYNoe8paLoawNT3hQFBkxTUHSUCFLEvrSBaiyhO7mww9fhbbvH8E3//VN5G0Xpy9owh2Xn+z3CNXYDunlNBgaOpoM7B7IY89gPhrGaoprGM4DW3vTkCChM1iVmSqv9s86IiKasuG8jZ0Hs0gaGtRgSCquK6NqamzXg2l7GMhacD2BloQGQ516zc3OgyP4xr++iazl4pSuFL5xxSmIaQosp/Z2SC8npinobopj90AeuwdzOA/zoo81xTUM5gS29vo9Pe0php5qULuDqERENC2O62HngRFYrodkrPw+VpoiozGmoq3RQEcqNq2w896hLP78p29ixHSwrCOJtVeeEoWpgu3W3A7p5cQ1BT3BTK09A/lxH29J6JAlCe/0prE/U5jr5tEEaj9mExHRlOwbyqMvXUBHcnZ6HHYP5PDnP30T6YKDE9obcecnTh3Vm1NwXHSmYjW3Q/pEdFXG8fPDmVq5CT+npUHHoRET7/SmIUsS2hpH1z5ZjgfTcWE6HkzHQ9504AqBRa0NNT2LrVox8BARHQPCoaxUrDiUNZP2DuZxx0/fwFDexpK2Btz1idPGrc1Tqzukl3NSZwpAsPigEBPOcpvXaEShZ8l8fx2iEdPBiOnAtF1YrgfbFQAEZMhwhIdMwcGJnUmkJumFo+lj4CEiqnOlQ1lTnWE1Hb3DftgZzNk4bl4Cf/HJ08atpOwJEeyQXj89Fx/oaIQsAXnbxaGsNa4HJxSGnrf2DkEAUCUZmipDV/zZcaosRWHJEwL96QLe2juMZZ0ptDboc/gT1TcGHiKiOjebQ1n96QK+/tibOJS10NOawF2fPA1NE/TihDuk10PBcigV0zA/aaA/bWLPYL5s4AH80DMVsiShMxXDgRETb+wdwkmdKXSw6HlG1M+ZR0RE4wznbew8NHNDWUIIvHcoh1ffH8Sr7w/irX1pOJ7AguY47v7kaWhOTNwjUcs7pJcT0xQsaI6jP21iy+4hfHBh04xsKyFJEtqTMQxmLby5bxiW43LLihnAwENEVKcc18Oug1lYjofW1JEPZaXzNv5rzxA2vzeI13YPYSBrjfr4kvkN+OYVp6BlkuEX0/bQ1RSrq4u2oco4fUETXn1/CP/31T3YeXAE//PiE2Zs7Z2WBh2Zgo13ejMwHQ/HtzXWRcF3pTDwEBHVqX1DefQO58cNZb3Tm8ZDv9yJEcuFrkjQlWJNSbjysq7IUGQJOw6O4Lf9IyjdB1wPLvRnLmrGmYtasKA5PmmQGSk4cIVX0zukT0SSJPyPcxZhKG/jiTd68er7Q1j9w1dx3TmL8KkzumekRy0Z06DIErbv92uwTpifrOltOSpJEtzPHgCQTqfR1NSE4eFhpFKpSjeHiOioDOdtbNk9CE2WR62588zWftz3zHY43vSe+he3JvA7i1pw5qJmnNrdNKWLbs5yMJizENcVLGxOYNG8xJT236ole4fyeGPvEOBJuP+57Xh97zAA4Lh5Caz+8AnRTK6jZTkeDoyY6GqKYVlnErE6WMtopkz1+s3AE2DgIaJ64bge3tqXRn+6EO3n5AmB//Pye/jJq3sAACuXzMPHl3fBdoU/NdrxYDkeLNcr3nY9tCcN/M6ilkkLcscq2C4GchYMRUZ3cxzdLXE01sF2EhMZyFr4z10D6EjFIAF4Zut+fP+XO5EpOJAAfOz0LnzuvMXjpuiXKtgueofz2DdUQEcqhhOCjUnHclwP/RkTrQ0a5jca0DUFmiLBUBRoqt9TNxtLDlQ7Bp5pYuAhonqQKdjoGyrg3YNZdCQNqIqMvOXifz25Df9v5wAA4Jqze/CZFYumvC/WVNmuvx0FJKC7KYYFLYkJZ2zVk7zl4tc7DwUF2X6vy3DexoO/3Ilngl3iWxt03HzhEhw3rwF7h/LYN5T33w777x8cGV0TdXJnEp88YwHOWzJvXM2OJwQGsxYs14OAACBBkSSoigRVkYOtOxQkdBWKLMHzBDwBCAj/fQCeV7ydimvobJreatrVhoFnmmoh8HiegCsEZEmCBEBm8RrVGcf14AoBIQA3ON89T0Tvh89WyZhaV9Obj5bnCQzlbfSl89ifNlGwPbQkNCR0FfszBfzlz9/BzoNZaIqEP/rIB3DJsvYZ/f6O62EgZ8ETQHvSQE9rAi0Jra4KlMsRQuD/7RyA44px4e6/9gzh/me3Y9/w4beWSMZUdKRi2HUwGw03ticNfHx5F/7bKZ2T9pC5noDtenA8ASdYyNDxPAgBSBIgQULwDxIQvC9BkgDT8dAc17BoXgLtQUCuNQw80zSbgad3OA/HFVBkCaosBW9lKIp/W5b8t54Qxe7l4L/leMhZLvK2A9MScCH8sBOcvbIEKLIEWQYUSYYkBfvgGP6rjbiuIKYeWTen5wn/j+UYeNKqB54nkLNd5EwHtucv8qZIEuTgHPPf99f5UGT/yS4MFp4Q8Dz/1aMniuFCCH/PoIShzFjthe36y+iHy+pbjodssPKs7YZtCdsBiODVqP9a1n8/pipoTuiYnzTQFNcmHS6oZ47r4VDWQm/QS+AJ/6IbhsGtvWnc/e/vYChnozmh4Y7LTz7impLoour6b23Piy7MsiShtUFDT2sCbQ3GMfdi7O3eYbx3MId5Dca4LSEsx8OP/3M3HnttLyQJ/hBfcxwLmuNY0BzzbzfFoxWoB7MWfv5mL/79jV6kCw4A/2/woye348rl3ehujs9o2z0hMJyzkXdctDXqWNTagHkN+owdQzdYWTpnOciZDha2Jma8N4mBZ5pmM/C8suMQDo6Y0BQpeMIuXnyU4OKjSDI8BOncE3Bdz4/mAFTJ76pUZf/ihbB7MjhynhCACN7CP8GcIKxoigRdVdBoqEjFVMR1f+M+XZXhitGvBhxXoGC7yNv+3i6240GWJcQ1BTFNRlxToSoSNEWGFnSfasEMj2MlFIW9bFFICH7vahBio2M0B0zHRd5ykbVcpPM2hnI2CrYDK1imPnhdFy15L0sSJBlQ4L+VIEGEvSkQEEHgCZ8Qwp9CliTEdQXNcQ3NDToadRUNhjJpiPY8P7gXgnPJtD2MmDayVnhuubDd4vdS5WL9gSz539P/jyiwhYQQyNsuRkwHjushpk0//Ph/I/4rYC86nv5jh8fUK+lRkoDgdzj6rRS0VULxdxcFNFFyX/A4ejgDSpWPanpxwXZxcMRf7G44Z0FVZDTH9VGFxM9t24///cxvYbsCx7c14M+vOBntwWwt1xNI523kbbd4oAEg6BGIjoyQoo+HwyaaIsNQZTQYKhK6Es3oak7ox+yU6eGcjV2HshjMWSjYHgy1+KIz5AbPyVMdRjQdF8//5gD+dcs+vD/g79UlATj3+FZ89KR2dDfHMT9pzFhPp+sJDGQtuEKgI2VgUWui7JpKh3ucrOUgazoYztkYzNrIOQ5sx4OqyDj3+NYZ3zKDgWeaZjvwjO3uDC+Yrlf8L5c8oczEE4cIe4zCV9KuByEEFNkPKv73L7nICf/iUuyFkuAJ/1WkE1zoQ35QQxDEZDQYCpKGCkNTYKhy9HaiXoGJul9t14M7ZtZI6fNCaaDygs8r/f05noDref7FygM8eFEPRfCjRRf30u+iyOHwYLGHLLzght/Tb6sH2/F7QVwUeyDC310YXFVJhqpK/u9AlRHTlGKPXvA7lUt+v2HvXmlIEqL0ZwoCllsMWnnLwWDORtZ0ULD9ISBFkhALgulEr54mupB7wfCoXHLhDi/kpRzXQ972w5XjCSiKH4JbEhqaEn4A8oRAwXFRsFxkTAdZ04XlurAdAVd4gJBKwnI4/Vk66qBcLvy0JY3oXLGDIlzLCXtOSwKrF/wuwhcQJS8cwh6l0vfDt+FwgB8cix8LHwNA9Hv2AwSA4G9vbGho0FXoqv/70FTZb0/pcF547ILzIQw7I6aDhKYiFddGPV9kTQcbX9uLH//nbgDAiuNb8dX/tgxxXYHpuBjO2XCFQHNCQ3syNuoCPNHhkGUpKowN21mLwx5zIWs6GMrb2J8uYChnwXQ9xBQFyZhWdlZb+PdefCGFUc+dQgj8155h/OuWvfjP9wbHfX2joaI9aWB+0kB70kB7Mob5SQPzGnSoigxFBhRZjp6jZBnR85EW/C2X/h1ajoeBnAlVlrCgOYEFLfFRLyKEKD43OZ7/3OR4/t9YOm9jYMRGPgg4siz7PcS6AgnAoZzFwDNT1q9fj7/6q79CX18fPvjBD+K+++7DueeeO6Wvna3A88QbvXh9zxAkSEjFtaAnoNhjoyrFiyFQvCCLkidNBLcFgpPM9U+0sFcmDAthgJBl/4Krq8qoV5R68CSryP4r+zBgKbJcfJEnFV/NSvADQelQXOl+L05wooeBxXa96IKghU/qmt+7pCkS8paLQtBzFAYUx/UvNr6SU1FIwc/tPxE4wc8ngqGXsJdLkiQIADLCsWq/7eEaIpoqR20Hxl/MRckrcFFyoYtuh8NCUti75hf9hT+v/6otKBoMnkyk4ArnBBfVseQgLIZDS+FbTZXgefAvyOGwjifgiLD3xfN/XgFowbEQQVvytouC7aJg+70qpaGm+DOVBD6B4ELqH8extTJRCJclJIKh0fCJS1dkSEFQ1GQZuibDjkK1iMJV+DYMHK4nRvWQOV7p9wxqdoJzuviE6gUBufg1flBQ0KCrfmAwVDQaKhp0f9hNChKIJzx4nn+uhyfmqACMYq1Q+LsKh/bCc84NEkz48xQ/Lww0/vthTV04hOj32iI6b2TJv7iEf+v+cS8OOWpBQNY1OUhO0Z9AFEIkAUiyBNf1X0FnCjb2ZyzsTxfQnylgf9pEf6aArOlG59rVZy7EZ89bhLztX4w0VcK8BgNdzTG0JnQGl1kihD+MM5Sz0Z8uYLhgw3YEdEWO/tYghSclovKG8DnWdFw4rgcJEgxViV487RvO4/HXe7G1N439GT/0Hi1DldHW6AemtkYdbY0G2hqNYGVuYEFzAp1Nhv+iwfZghdcf1w3+Rka/gFcVGVrwfBs+bzvBC45DOQufPW8ROlMzOyx3zAWeH/3oR/jc5z6HBx54ACtWrMC9996LRx99FNu2bUN7++EL9GYr8Hzkfz2HHQeyM/Z4lTY2AIUX7LB3IOwxCbv7w/fHnmZRxAkvQGGocYVfG+AWe4FmghzUNoWvrlVZLunJwYS9HBKK9SZmSb3JVJqkylIUMsPQVRpwwqGaMFSGhYVRj58YHQ5KA20hqH8hKmdeg47PnrcYZy1uQc5ykdAVdDXF0J6MIRVXj5kh6GrgeQIZ08FQzkKmYAcBRgle7EpBL4wfiMPnh7BH1e8xspA1XZhBb64s+S9oVdkPRgNZG4eyJg6OWDiYMXFgxMT+tInBnDVuFKH0+WSGnlqn7fEvfwinLWia0cc85gLPihUrcM455+Dv//7vAQCe56Gnpwdf/vKXcdttt437fNM0YZpmdDudTqOnp2fGA8/XH3sD7+xL+6+6EbyKDYeIvGLPRTicE4YEoHgRBopV9uqYwKEGXeSKJEWvElzPH8ayPRGtpWGXdOlbTjA8JYrFoBgz1FPNirMMxsw+CH5XIuhdmW2y5K84C1HslZlrkoTo1V9MVWBoMgxF9sNUyUyM0nMq/LpiAMOonqbSt64QMIO6rrD3qDDmtgh/FyWr9ZaGvXD4KhrOk8a+LdbplAbp8H05HBIM2mW5fiF/znKit3nLRc5ykQ3us12v+LNEvSwY9b2l4O8p7JEJf+7Rv4PRf4cySv4mg1AcjmeFPT8irIfC6J6ucPjWDOqXwmG28O9yOudso6FEr8TbGnXMa/SHMFobDbQldGiqDECgKaFjQXMcbY3ji2mpdoT1erkoBNmwHa+k/mx0eUI4xOpXjPq9wqXjs5JU7KUfztsYzFkYyFoYzFoYyNnBW/++nOWWadV4YWgL6xnDEQxV9q9PsiRh/Wd+B6d0VSbw1MXUBsuysHnzZtx+++3RfbIsY9WqVXj55Zcn/Jp169bhzjvvnPW2fft3T8e2vjQyQbV96dBCeBso9nRExNibxTvCizxQOoxT8gQ8geL3meR7R/eLYCgJZV8dFIcbwif00id3f1igdJggbJv/ylKMnioZvB/+UeiyXOyNUYuF0eEw4OEI+IHXcQWs4MLiOMVxZjsoDA+HrUrbOeo+iKDWQolqcvz/CnTNb2Npa5ygXiQs+LYd4dexuH438Nghk7CGxA3u8yBGz+KTxvemqbIUhZuY7i86Vlr/U3ouhLejj0kYHaBRDEQY83mlj1d6KhZ75Py3nvDrm1RF9h+r5DuOrsEa/xgTHjtR0saorRgXMsK/o/BrJmrjVIz9WceFmeAjpW0Ka3aizysJ2tEQaXQ7atWY4enRbfVDugfbc/3hRq/kMaKaMf+2GtRdhOeJLBcDa+m5EtcUtDTodbey8bEo7BVqThTv88rU+TmeB89DVJQfCs+l8P3o+Q7F5z0Po89jz/Nr5By3+JwdPY9L/pkeDrMXn9/HX2/C27oi4/h5Ey+qOBfqIvAcPHgQruuio6Nj1P0dHR3YunXrhF9z++23Y82aNdHtsIdnNiyb4jTQiTrbyj15l17UZkLpxQOYuLfncJ2Bpe0pXhDGPkbJ+xM8djEYjX/MozE26JXW64z/3NG3S4fmptK2yX6X0+lQnejxx/5eOTRBdGySZQn6FF4AUlFdBJ4jYRgGDOPIdw+eDRNe4ObofI56BCb9fkffmPKPP7s/aOkr8dn+fpP/LvkERURUCXXR19nW1gZFUdDf3z/q/v7+fnR2dlaoVURERFQt6iLw6LqOs846C08//XR0n+d5ePrpp7Fy5coKtoyIiIiqQd0Maa1ZswY33ngjzj77bJx77rm49957kc1mcdNNN1W6aURERFRhdRN4rr32Whw4cADf/OY30dfXhzPOOAObNm0aV8hMREREx566WYfnaNXCbulEREQ02lSv33VRw0NEREQ0GQYeIiIiqnsMPERERFT3GHiIiIio7jHwEBERUd1j4CEiIqK6x8BDREREdY+Bh4iIiOpe3ay0fLTC9RfT6XSFW0JERERTFV63D7eOMgNPIJPJAAB6enoq3BIiIiKarkwmg6amprIf59YSAc/zsG/fPiSTSUiSVOnmVJ10Oo2enh7s3r2bW29UCR6T6sLjUV14PKrLbB4PIQQymQy6u7shy+UrddjDE5BlGQsXLqx0M6peKpXik0eV4TGpLjwe1YXHo7rM1vGYrGcnxKJlIiIiqnsMPERERFT3GHhoSgzDwNq1a2EYRqWbQgEek+rC41FdeDyqSzUcDxYtExERUd1jDw8RERHVPQYeIiIiqnsMPERERFT3GHiIiIio7jHwUGT9+vU47rjjEIvFsGLFCvz6178u+7nf+973cOGFF6KlpQUtLS1YtWrVpJ9P0zed41Fqw4YNkCQJn/rUp2a3gceg6R6ToaEhrF69Gl1dXTAMAyeeeCKeeOKJOWpt/Zvu8bj33nuxbNkyxONx9PT04NZbb0WhUJij1ta3F154AVdeeSW6u7shSRJ++tOfHvZrnnvuOZx55pkwDAMnnHACHn744dltpCASQmzYsEHoui4efPBB8dZbb4kvfOELorm5WfT390/4+ddff71Yv369eO2118Q777wjfv/3f180NTWJPXv2zHHL69N0j0do586dYsGCBeLCCy8Un/zkJ+emsceI6R4T0zTF2WefLS6//HLx0ksviZ07d4rnnntObNmyZY5bXp+mezweeeQRYRiGeOSRR8TOnTvFf/zHf4iuri5x6623znHL69MTTzwh7rjjDrFx40YBQDz22GOTfv6OHTtEIpEQa9asEW+//ba47777hKIoYtOmTbPWRgYeEkIIce6554rVq1dHt13XFd3d3WLdunVT+nrHcUQymRT//M//PFtNPKYcyfFwHEecf/754p/+6Z/EjTfeyMAzw6Z7TL773e+KJUuWCMuy5qqJx5TpHo/Vq1eLj3zkI6PuW7NmjbjgggtmtZ3HoqkEnj/7sz8Tp5566qj7rr32WnHZZZfNWrs4pEWwLAubN2/GqlWrovtkWcaqVavw8ssvT+kxcrkcbNtGa2vrbDXzmHGkx+Mv/uIv0N7ejs9//vNz0cxjypEck5/97GdYuXIlVq9ejY6ODpx22mn49re/Ddd156rZdetIjsf555+PzZs3R8NeO3bswBNPPIHLL798TtpMo7388sujjh8AXHbZZVO+5hwJbh5KOHjwIFzXRUdHx6j7Ozo6sHXr1ik9xte+9jV0d3ePO4Fp+o7keLz00kv4/ve/jy1btsxBC489R3JMduzYgWeeeQaf+cxn8MQTT2D79u340pe+BNu2sXbt2rlodt06kuNx/fXX4+DBg/jQhz4EIQQcx8EXv/hFfP3rX5+LJtMYfX19Ex6/dDqNfD6PeDw+49+TPTx01O655x5s2LABjz32GGKxWKWbc8zJZDK44YYb8L3vfQ9tbW2Vbg4FPM9De3s7/vEf/xFnnXUWrr32Wtxxxx144IEHKt20Y9Jzzz2Hb3/727j//vvx6quvYuPGjfj5z3+Ou+66q9JNoznCHh5CW1sbFEVBf3//qPv7+/vR2dk56df+9V//Ne655x489dRTWL58+Ww285gx3ePx7rvvYteuXbjyyiuj+zzPAwCoqopt27Zh6dKls9voOnckfyNdXV3QNA2KokT3nXzyyejr64NlWdB1fVbbXM+O5Hh84xvfwA033IA//MM/BACcfvrpyGazuPnmm3HHHXdAlvn6fy51dnZOePxSqdSs9O4A7OEhALqu46yzzsLTTz8d3ed5Hp5++mmsXLmy7Nd95zvfwV133YVNmzbh7LPPnoumHhOmezxOOukkvPHGG9iyZUv0/xOf+AQ+/OEPY8uWLejp6ZnL5telI/kbueCCC7B9+/YofALAb37zG3R1dTHsHKUjOR65XG5cqAnDqOCWknNu5cqVo44fADz55JOTXnOO2qyVQ1NN2bBhgzAMQzz88MPi7bffFjfffLNobm4WfX19QgghbrjhBnHbbbdFn3/PPfcIXdfFT37yE9Hb2xv9z2QylfoR6sp0j8dYnKU186Z7TN5//32RTCbFLbfcIrZt2yYef/xx0d7eLv7yL/+yUj9CXZnu8Vi7dq1IJpPihz/8odixY4f4xS9+IZYuXSquueaaSv0IdSWTyYjXXntNvPbaawKA+Ju/+Rvx2muviffee08IIcRtt90mbrjhhujzw2npf/qnfyreeecdsX79ek5Lp7lz3333iUWLFgld18W5554rXnnllehjF198sbjxxhuj24sXLxYAxv1fu3bt3De8Tk3neIzFwDM7pntMfvWrX4kVK1YIwzDEkiVLxN133y0cx5njVtev6RwP27bFt771LbF06VIRi8VET0+P+NKXviQGBwfnvuF16Nlnn53wmhAegxtvvFFcfPHF477mjDPOELquiyVLloiHHnpoVtsoCcG+PCIiIqpvrOEhIiKiusfAQ0RERHWPgYeIiIjqHgMPERER1T0GHiIiIqp7DDxERERU9xh4iIiIqO4x8BAREVHdY+AhIiKiusfAQ0RERHWPgYeIiIjqHgMPEdWlH/7wh4jH4+jt7Y3uu+mmm7B8+XIMDw9XsGVEVAncPJSI6pIQAmeccQYuuugi3HfffVi7di0efPBBvPLKK1iwYEGlm0dEc0ytdAOIiGaDJEm4++67cfXVV6OzsxP33XcfXnzxRYYdomMUe3iIqK6deeaZeOutt/CLX/wCF198caWbQ0QVwhoeIqpbmzZtwtatW+G6Ljo6OirdHCKqIPbwEFFdevXVV3HJJZfgH/7hH/Dwww8jlUrh0UcfrXSziKhCWMNDRHVn165duOKKK/D1r38d1113HZYsWYKVK1fi1VdfxZlnnlnp5hFRBbCHh4jqysDAAM4//3xccskleOCBB6L7r7jiCriui02bNlWwdURUKQw8REREVPdYtExERER1j4GHiIiI6h4DDxEREdU9Bh4iIiKqeww8REREVPcYeIiIiKjuMfAQERFR3WPgISIiorrHwENERER1j4GHiIiI6h4DDxEREdW9/w/XDfGpAYIVSgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Background\n",
"plt.plot(x, samples_bkg.mean(0), color=\"C1\", label=\"Background\")\n",
"plt.fill_between(x, \n",
" np.percentile(samples_bkg, 16, axis=0), \n",
" np.percentile(samples_bkg, 84, axis=0), \n",
" color=\"C1\", alpha=0.2)\n",
"\n",
"# GP\n",
"plt.plot(x, np.exp(samples['log_rate'].mean(0)), color=\"C0\", label=\"GP\")\n",
"plt.fill_between(x, \n",
" np.percentile(np.exp(samples['log_rate']), 16, axis=0), \n",
" np.percentile(np.exp(samples['log_rate']), 84, axis=0), \n",
" color=\"C0\", alpha=0.2)\n",
"\n",
"# Combined\n",
"plt.plot(x, (np.exp(samples['log_rate']) + samples_bkg).mean(0) , color=\"C2\", label=\"Combined\")\n",
"plt.fill_between(x, \n",
" np.percentile((np.exp(samples['log_rate']) + samples_bkg), 16, axis=0), \n",
" np.percentile((np.exp(samples['log_rate']) + samples_bkg), 84, axis=0), \n",
" color=\"C2\", alpha=0.2)\n",
"\n",
"\n",
"plt.errorbar(x, y, yerr=y_err, fmt='o', color='k', label=\"Realized counts\")\n",
"\n",
"plt.xlabel(\"$x$\")\n",
"plt.ylabel(\"$y$\")\n",
"\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16206320-65b0-43d2-ac6f-372f711b36ff",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "bfa81cff-d106-4981-ab12-f8c76bdd01a1",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment