{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multibin Coupled HistoSys\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import logging\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pyhf\n",
    "from pyhf import Model\n",
    "from pyhf.contrib.viz import brazil\n",
    "\n",
    "logging.basicConfig(level=logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep_data(sourcedata):\n",
    "    spec = {\n",
    "        \"channels\": [\n",
    "            {\n",
    "                \"name\": \"signal\",\n",
    "                \"samples\": [\n",
    "                    {\n",
    "                        \"name\": \"signal\",\n",
    "                        \"data\": sourcedata[\"signal\"][\"bindata\"][\"sig\"],\n",
    "                        \"modifiers\": [\n",
    "                            {\"name\": \"mu\", \"type\": \"normfactor\", \"data\": None}\n",
    "                        ],\n",
    "                    },\n",
    "                    {\n",
    "                        \"name\": \"bkg1\",\n",
    "                        \"data\": sourcedata[\"signal\"][\"bindata\"][\"bkg1\"],\n",
    "                        \"modifiers\": [\n",
    "                            {\n",
    "                                \"name\": \"coupled_histosys\",\n",
    "                                \"type\": \"histosys\",\n",
    "                                \"data\": {\n",
    "                                    \"lo_data\": sourcedata[\"signal\"][\"bindata\"][\n",
    "                                        \"bkg1_dn\"\n",
    "                                    ],\n",
    "                                    \"hi_data\": sourcedata[\"signal\"][\"bindata\"][\n",
    "                                        \"bkg1_up\"\n",
    "                                    ],\n",
    "                                },\n",
    "                            }\n",
    "                        ],\n",
    "                    },\n",
    "                    {\n",
    "                        \"name\": \"bkg2\",\n",
    "                        \"data\": sourcedata[\"signal\"][\"bindata\"][\"bkg2\"],\n",
    "                        \"modifiers\": [\n",
    "                            {\n",
    "                                \"name\": \"coupled_histosys\",\n",
    "                                \"type\": \"histosys\",\n",
    "                                \"data\": {\n",
    "                                    \"lo_data\": sourcedata[\"signal\"][\"bindata\"][\n",
    "                                        \"bkg2_dn\"\n",
    "                                    ],\n",
    "                                    \"hi_data\": sourcedata[\"signal\"][\"bindata\"][\n",
    "                                        \"bkg2_up\"\n",
    "                                    ],\n",
    "                                },\n",
    "                            }\n",
    "                        ],\n",
    "                    },\n",
    "                ],\n",
    "            },\n",
    "            {\n",
    "                \"name\": \"control\",\n",
    "                \"samples\": [\n",
    "                    {\n",
    "                        \"name\": \"background\",\n",
    "                        \"data\": sourcedata[\"control\"][\"bindata\"][\"bkg1\"],\n",
    "                        \"modifiers\": [\n",
    "                            {\n",
    "                                \"name\": \"coupled_histosys\",\n",
    "                                \"type\": \"histosys\",\n",
    "                                \"data\": {\n",
    "                                    \"lo_data\": sourcedata[\"control\"][\"bindata\"][\n",
    "                                        \"bkg1_dn\"\n",
    "                                    ],\n",
    "                                    \"hi_data\": sourcedata[\"control\"][\"bindata\"][\n",
    "                                        \"bkg1_up\"\n",
    "                                    ],\n",
    "                                },\n",
    "                            }\n",
    "                        ],\n",
    "                    }\n",
    "                ],\n",
    "            },\n",
    "        ]\n",
    "    }\n",
    "    pdf = Model(spec, poi_name=\"mu\")\n",
    "    data = []\n",
    "    for channel in pdf.config.channels:\n",
    "        data += sourcedata[channel][\"bindata\"][\"data\"]\n",
    "    data = data + pdf.config.auxdata\n",
    "    return data, pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "validation_datadir = \"../../../validation/data\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pyhf.pdf:Validating spec against schema: model.json\n",
      "INFO:pyhf.pdf:adding modifier coupled_histosys (1 new nuisance parameters)\n",
      "INFO:pyhf.pdf:adding modifier mu (1 new nuisance parameters)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data: [110.0, 105.0, 170.0, 220.0, 0.0]\n",
      "parameters post unconstrained fit: [-0.30257894  0.63607078]\n",
      "parameters post constrained fit: [0.29087602 0.        ]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([106.80676913, 104.01075137, 154.36314037, 213.86871072,\n",
       "         0.29087602])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open(\n",
    "    validation_datadir + \"/2bin_2channel_coupledhisto.json\", encoding=\"utf-8\"\n",
    ") as spec:\n",
    "    source = json.load(spec)\n",
    "\n",
    "data, pdf = prep_data(source[\"channels\"])\n",
    "\n",
    "print(f\"data: {data}\")\n",
    "\n",
    "init_pars = pdf.config.suggested_init()\n",
    "par_bounds = pdf.config.suggested_bounds()\n",
    "\n",
    "unconpars = pyhf.infer.mle.fit(data, pdf, init_pars, par_bounds)\n",
    "print(f\"parameters post unconstrained fit: {unconpars}\")\n",
    "\n",
    "conpars = pyhf.infer.mle.fixed_poi_fit(0.0, data, pdf, init_pars, par_bounds)\n",
    "print(f\"parameters post constrained fit: {conpars}\")\n",
    "\n",
    "pdf.expected_data(conpars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert_interval(test_mus, cls_obs, cls_exp, test_size=0.05):\n",
    "    crossing_test_stats = {\"exp\": [], \"obs\": None}\n",
    "    for cls_exp_sigma in cls_exp:\n",
    "        crossing_test_stats[\"exp\"].append(\n",
    "            np.interp(\n",
    "                test_size, list(reversed(cls_exp_sigma)), list(reversed(test_mus))\n",
    "            )\n",
    "        )\n",
    "    crossing_test_stats[\"obs\"] = np.interp(\n",
    "        test_size, list(reversed(cls_obs)), list(reversed(test_mus))\n",
    "    )\n",
    "    return crossing_test_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_pois = np.linspace(0, 5, 61)\n",
    "results = [\n",
    "    pyhf.infer.hypotest(\n",
    "        test_poi, data, pdf, init_pars, par_bounds, return_expected_set=True\n",
    "    )\n",
    "    for test_poi in test_pois\n",
    "]\n",
    "cls_obs = np.array([result[0] for result in results]).flatten()\n",
    "cls_exp = [np.array([result[1][i] for result in results]).flatten() for i in range(5)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'exp': [0.3379557904310767,\n",
       "  0.6355433028074602,\n",
       "  0.980748685945719,\n",
       "  1.4196566943996904,\n",
       "  1.9446272307776151],\n",
       " 'obs': 1.5452056912523955}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFFCAYAAAB47zMwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABbX0lEQVR4nO3dd3hURffA8e+kF0KAUAKEXhOqEKqAgNSI+EMFAUUpFhBUwEJRUFBBVBRUQBGV8irVAoIC0gQF6VV6FYJ06Wkk8/tjNrCE9OzdTcj5PE+e7O6de+ewry+HuXdmjtJaI4QQQuQGbq4OQAghhHAWSXpCCCFyDUl6Qgghcg1JekIIIXINSXpCCCFyDUl6Qgghcg2nJT2l1NdKqTNKqV0pHFdKqU+UUgeVUjuUUrWcFZsQQojcwZkjvalAm1SOtwUq2H6eBSY5ISYhhBC5iNOSntZ6NXAhlSYPAdO18ReQTylV1DnRCSGEyA2y0zO94sBxu/cnbJ8JIYQQDuHh6gAyQyn1LOYWKP7+/rUrV66c6WtduXKF/fv3A1CjRg08PHLkVyKEEMJm8+bN57TWhZI7lp3+ho8ESti9D7F9dget9WRgMkB4eLjetGlTpju9ePEi+fPnp0JF2LatF/ACWmuUUpm+phBCCNdRSh1L6Vh2ur25AHjSNouzPnBJa/2v1Z3my5ePkFIFKVYB4EUWLXqABx9sx+XLl63uWgghhJM5baSnlJoJNAUKKqVOAG8CngBa68+BX4AI4CBwHejhrNiKBBdk975zHD0KZ8/+wrlzgbi7XwHyOisEIYQQTqByemmhrN7eBG7eynz+eZgwAeLjwd29MgkJ8zl9OoCiRWUSqRBC5BRKqc1a6/DkjmWnZ3ou0/De2qw9sJluz5j37u4Ae3nmmVqsXl2Ybdt24u/v78oQhRDZXFxcHCdOnCA6OtrVoeQaPj4+hISE4Onpme5zJOkBf/6xBP9RBTnhdfvnPXteIzzcGz8/H9cEJoTIMU6cOEFAQAClS5eWiXBOoLXm/PnznDhxgjJlyqT7vOw0kcVlzp+PwScSpnwHmzff+vzee6FPn70oNYzr16+7LkAhRLYXHR1NUFCQJDwnUUoRFBSU4ZG1JD0gNLQmFybCb2Ph00/vPL5z52jKlw9hyZIlzg9OCJFjSMJzrsx835L0gCeffBKAhDowZNSdx8uWhXvvvULJkglOjkwIIYQjSdIDxowZg3IDFKw7c+dxf3+YO/cGoaGvAtfI6TNehRAit5Kkh9mKzCsA+A8mfAaLF6fU8m+GDq3D008/LYlPCJEtnTp1is6dO1OuXDlq165NREQE+/fvJ0+ePK4OLVuQpIe5vRlzCbgAW3+An35KuW2+fHvw9DwkO7YIIbIdrTUdOnSgadOmHDp0iM2bNzN69GhOnz7t6tCyDVmyAHTv3p1tR5dwvFws8ZXgpZdTbvvaawD+QKCTohNCiPRZuXIlnp6e9O7d++ZnNWrUcGFE2Y+M9ICHH36Ylz4vBzUBX1h+JK0z/uD8+WQe/gkhhAvt2rWL2rVruzqMbE1Gephneie23IADQCn4aBjk6QLduyff/q23LvP++yW5ePEyXl5eyTcSQuRa/fv3Z9u2bQ69Zs2aNRk3bpxDr5kbSdIDZsyYwbgXDpg3feGfPXDgYMrtW7eGwMBWxMXFSdITQmQbVapUYd68eeluf/z4cUaOHElgYCBt2rShRYsWFkaXPUjSA1q0aEHd+4qw4ffTcAXie8L/PZNy+wYNoEGDBMyzPSGEuJ2rRmTNmzdn6NChTJ48mWeffRaAHTt2cOnSpWTb7927Fy8vL1588UVKlizpzFBdRp7pARUrVmToqMbmzWVApf1cLypqDVu3Zq26gxBCOJJSih9//JFly5ZRrlw5qlSpwpAhQwgODub69euEhITc/Pnoo49o2bIlL7zwAv369SMyMtma3XcdGekBMTExnD9u223lMpAAHz4P8Y/B668nf87gwZeZMqUxFy9eztAO30IIYaVixYoxZ86cOz5PSLhzR6lBgwYRHx9PyZIlKVy4sDPCczlJesD+/fvp1fkHUJik5wYX3cE/lbWcTz8Nbdo87qwQhRDC4caMGePqEJxObm8CpUqV4tXX2kBDwDbbN/7/oEb7lM+pVg3atj0tozwhhMhBJOkBefPm5f0xPQhsB9gVSV9+GOLiUj7vwIFVLFjwk9XhCSGEcBBJepite1asOIT7dmCv7cN4GPsYjBiR8nmffHKVrl27EpdaZhRCCJFtSNKzuf/+oVyYA8wCbgDuEF8VatVK+ZyBA2H79v54eMijUSGEyAkk6WGm+X700UCK1Ld9cMX8imsGZRumfF6ZMlCu3C4pHCmEEDmEJD2bAQOeo+x9tjd2BRRW7oUrV1I+b8WKFUyaNMHS2IQQQjiGJD2bdev2EnPC9iYxyV2Egc3hu+9SPu+HH67x1lvDk10DI4QQInuRpGfz+OMvsmWG7U3iSC8Qgh6Ahqnc4hw5Eo4fH4ybm3yVQgiR3cnf1Dbjxn1MxEtAZ6Ca7UMF58OhaLmUzytQALy8/nBChEIIkTapnJ46SXo27dt3oHUXD6gMBNgdiIdvl8LVqymf+803yxg27A2rQxRCiFRJ5fS0SdKz2b17N5t/coNdmJ9Ex6F/B/j995TP3bLlOsuWLUBrbXWYQgiRopQqpzdu3NiFUWUvssDM5p133mHmzFizI8t1oKrtQDEI6w1166Z87rhx4O7eE7N5pxBCuIZUTk+bJD2bIUOGULjSYsYv+g82AwmYcbAX/FMSCgSlfK67O8AqoL/1gQohcoSmTZve8VmnTp14/vnnuX79OhEREXcc7969O927d+fcuXM8+uijtx1btWqVRZHmLpL0bKpVq0avwiGMX/+fSXjXAdtz36sX4MPJ8MoziQnuTm++uYSTJ5/myy+nOCtkIYS4jSMrpx87doyxY8eitaZcuXL079+fH3/8kePHjxMWFsaGDRsoVqwYK1euJCwsDHd3d1555RUr/lgOJUnP5vjx48ydcQVibR9c5mbS4xAM/hAi7jXVFZJz40Y0cXFnnBCpECInSG1k5ufnl+rxggULZmpk58jK6RMnTsTX1xdfX1927twJQIcOHejXrx9r1qxh1qxZzJgxg5YtW/LEE0/QuXNnEhISsv3yLUl6NvPmzePt149CcdsHl4FittcV4IEPITQ05fPffRfAF4ghPt4D95SGhEIIYZHEyun9+/dnzJgx+Pj4ULp0acaNG3ezcnqigQMHMnDgQEqUKEG/fv2YNGkSxYsXv3k8ISGBbt26Ub169Zufaa25ePEibm5u3LhxA+Dm75yy8b4kPZuOHTtSqNBcnly3Du0G5LM76AcHvCHtfaXnsH37bjp3jmbOnB+oltKwUAghLJKVyumnT59m4cKF9OrVi379+jF06FCKFi1KQEAAb775JuPHj6dLly4EBwczbNgwwsLCWLp0KTt27KBu3brZfpQHoHL6NPvw8HC9adMmB13tOQp9MJlz15M59C+MrgCDX079Ctu2Qd++nsyb9wNFi7ZzUFxCiOxuz549hKZ2OygHWLx4MV5eXjRv3jxd7adOnUrBggVp1851f9cl970rpTZrrcOTay8jPZtz584xa9YB3LYDVzHfjP1A7SC8/iX0eRoCA1O+Ts2a8McfcSj1KFpPZNYsXzp16iS3O4UQ2V6bNm0y1L579+7WBGIhSXo2O3fu5IUXVqI8gGDAk9uTXjgMeDH1hJfIVBqKYfnyXnTtCgkJcTz++JNWhC2EECIDsv8NWCepXbs2ixb1oO5wIC+3lRcCwBe2XszYNVu0gCVLoGvXHzGVaYUQQriSJD2bvHnzEhFRh5KVuZX0kjzuXLsYRo3K2HVbtQKlfuLkycd47bVXb850EkII4XyS9Gyio6P58suNHFsE+ABxQEySNgdgylTIzNyfxYt/YNKk8ezbty/rwQohhMgUSXo2169f59lnv2HDNG5toXkRiLdr1AZe+zbxmV3G9OwJ+/bFUaXK1KyGKoQQIpMk6dkEBgaydu1H9J0C1AOeAqYAx+waecBfkZnvo1gxgA+ZPbsTr7/+ulRlEEIIJ5OkZ+Pu7k6DBvWoUh1ze7MwZu7Jv7e3+/kLeOedrPW1Zs1c1qyZRWxsbNqNhRBCOIwsWbAzd+5aFi8B3DE19fy4I+ldOAF79metn08+gZiYw3h7/wB0ydrFhBBCpJuM9OwMGPAhC74CDgHbAC/uSHo8Cl2zWCTdzQ18feHKlad57rnHOHr0aNYuKIQQIl0k6dlZteonpv4K1McsW3AHznPHLM5tpxzT39mz15k373tWr/7NMRcUQog09OzZk8KFC1O1atW0Gyfj1KlTdO7cmXLlylG7dm0iIiLYv9/c/sqTJ08aZ7ue3N60U758NW7kA9Zjkl5ipaBTQCnb68sw6UWoNRrats1af2XLwqFD8eTLtxZ4JmsXE0JkK2pEJqZ5p0K/mbGJb6tWrWLq1KlMnTr1ts+7d+9Ov379ePLJjO8SpbWmQ4cOPPXUU8yaNQuA7du3c/r0aSpWrJjh67mCU0d6Sqk2Sql9SqmDSqnByRwvqZRaqZTaqpTaoZS6s7SwhZYtW8snbwN7gABujfDsb3F6wfXr4KgqGvnyAUxl06Y3WbRokWMuKoQQKWjSpAkFChS44/Pt27fTpEkTwsLCcHNzQynF8OHDb2uzcuVKPD096d27983PatSoQePGjS2P21GcNtJTSrkDE4CWwAlgo1JqgdZ6t12zN4A5WutJSqkw4BegtLNiHD/+E1aswOy9mTiy8+f2pOcD5V6F9u0d16/WMHDgO1y8WI62bdvmiPIcQojsqV69esTExHD16lUuXLhAzZo1ARgzZgytW7dO9pzo6Ggee+wxpk+fTt26dRk2bBjR0dGMGDHitna7du2idu3aVv8RLOXM25t1gYNa68MASqlZwEOAfdLTmBuLAIHASSfGx5dffsmlS1WpMfO8GeTdB3zLHZNZTl5xbL9KwXffJeDrC25u17lVsl0IITJm/fr1QMq3N5OzbNkyatWqRd26dQGoXr06ixcvRmVmJ45szplJrzhw3O79CcwycHtvAUuVUi9gxlgtkruQUupZ4FngjhL3WREcHExwcD4KBJzn36u2D4sCB4FYzGxO4N+50G0bzJjusK4xBY0PkJDwHMuXd6dly5aOu7gQQqRi165dtxW93rJlC7Vq1bqjXZUqVZg3b166r3v8+HFGjhxJYGAgbdq0oUWL2/9KP3bsGGPHjkVrTbly5ShVqhTHjx8nLCyMDRs2UKxYMVauXElYWBju7u688sormf9D2mS3+2hdgKla6xAgApihlLojRq31ZK11uNY6vFChQg7rfP369QwceIn4P4BrwBzMlmSaW5NaAO0NPhYNxr766jtatWqF4wrjCiFyo6ZNm6ZrlAcQFBTEjh07ANi/fz8//PADnTt3vqNd8+bNiYmJYfLkyTc/27FjB2vWrEn2unv37sXLy4sXX3zxjoQHMHHiRHx9fQkKCmLnzp106NCB/fv388UXXzBo0CAAWrZsyaBBg9i0aVOy1d8zypkjvUighN37ENtn9noBbQC01uuUUj5AQW5LOdb5+eef+fjjcyh34DXMjdd7bQf/tUUM0Ax6WzTZ8vHHwd8/L7Vrl7WmAyHEXS3xmV5Sic/0unTpwqpVqzh37hwhISGMGDGCLl26sGDBAqpWrUrBggWZOXMmQUFBd1xDKcWPP/5I//79GTNmDD4+PpQuXZpx48YBZg/jkJCQm+0HDhzICy+8QL9+/Zg0aRLFixe/7XoJCQl069aN6tWrA2Z26MWLF3Fzc7tZkSbxd5yDZg8qZ+3/qJTyAPYD92OS3Uagq9b6b7s2vwKztdZTlVKhwHKguE4lyPDwcO2oUdGlS5e4fLk9z6xYzZKjwBggDJP8QgG7ySsLOsODlRzSbQqe5saNSXh4yKoSIXKCPXv2EBoa6uowspVBgwYRHx9PdHQ0H3/8MRcuXGDhwoX06tULMLc3hw4dStGiRQkICCAwMJAKFSoQHBzM7NmzCQsLY+nSpQQHB1OkSJGboz97yX3vSqnNWuvw5GJy2t+oWusbSql+QOJGX19rrf9WSo0ENmmtFwAvA18qpQZgbip2Ty3hOVpgYCCBgYUICkj8AFNXryi3T2Y5DM+0gt+XQiWLEt/atVN4/PFFLFmyKsesfxFCCHtjxoy57f3WrVspU6bMzfelSpXi22+/Tfbc2rVrM3XqVLp27Uq7du0cFpNThxFa618wyxDsPxtu99r+hqLT7du3j6lT97H7AGbJQiDwH1ABWIfZgNoD8IWgspkrMZReZctC2bKXiI29Zl0nQgjhRG3atMlQ++7duzs8Brl3ZmfHjh28994u86YnptJCFGakl4B5sljMvG8QAVYOwIKDYfny68Ai4B7rOhJCiFwku83edKn27dtz/vwLjF4FlMQ8feyJSXRw2y1OR6/VS0lU1DuMGPES586dc06HQghxF5OkZ8fb25sCBQpRMCDJgfyAN7clvVVDYfAdG6k53uHDMbz77qfMnz/f+s6EEOIuJ7c37Zw9e5bPP1/L6h2YpfF5gR+B5tw5maUkOGOiVpUqsH+/pnRpb+s7E0KIu5yM9OxcvHiR4cMXs/wHzKIKT8y+MecxSe80EG/aRreArk84J67SpQEGcvLknptrVoQQQmScJD075cqV4/r1b1h/DGiE2QLTDbiESXo3ANujNQ23tipzgr//Pku5ctWZPt2Be58JIUQuI0nPjpubG76+BQjyT/wAc4szMenBrVucf0DVUqZCgjOEhcFrr0GzZrL4VQghMkue6dnRWvPOOz9y5DimvkMdzFq9S0AQ5nbnv0BNoAg0/T+IjQVvJzxuUwpGjLgBfAM0sL5DIYS4C8lIz45SitGjZ7JiKXDY9mEpTMJzw9TZSxzpVYBWzzsn4dk7ceIr+vTpwsWLF53bsRAixzt+/DjNmjUjLCyMKlWqMH78+Axf49SpU3Tu3Jly5cpRu3ZtIiIi2L9/PwB58mT/smiS9JK4fHkDR49CQDfbB80xVf/g1gxO20bfkZchPt658Z07l8C0aXNZt26dczsWQmSQcvBPxqxateqOHU08PDwYO3Ysu3fv5q+//mLChAns3r07+QskQ2tNhw4daNq0KYcOHWLz5s2MHj2a06dPZzg+V5Gkl4SHRyAA+X2TOVgUiAMuAJfggwdg2jQnBgfUrAmRkfG0bXvnDuhCCJGaokWL3qyTFxAQQGhoKJGRptjN9u3badKkCWFhYbi5uaGUYvjw4bedv3LlSjw9Pendu/fNz2rUqEHjxo2d94fIInmml8SECbPZsweu7cSM8o4Dc4FHuX0ySygUv99MMHG2/PkBBhMZOeOOUh1CiNwtsbTQ1atXuXDhAjVr1gRulRZKdPToUbZu3Uq9evWIjo7mscceY/r06dStW5dhw4YRHR3NiBEjbrv2rl27qF27tjP/OA4nI70kvv32R1asgKu7MLcxvTGVFi4BhTD1If4FPCBPO6hf3zVxTp26kjJlSt+8ly6EEGCKYW/bto0pU6bQvn17tm3bxrZt225LeFevXuWRRx5h3Lhx5M2bl2XLllGrVi3q1q0LQPXq1blw4QLKyl31XURGekn8+edalPLm/2bFMX8fZvYmmKTnDhTh5mSWyMtw/Tr4+Tk/zrZtYfDgIIoUKez8zoUQOVZcXByPPPIIjz/+OA8//DBgRnDVqlW72WbLli03b4Paq1KlCvPmzUt3X8ePH2fkyJEEBgbSpk2bO6qnHzt2jLFjx6K1ply5cpQqVYrjx48TFhbGhg0bKFasGCtXriQsLAx3d3deeeWVTP6pb5GRXhLmXzb+5POxfeAN+GKSHtyazKLh0jfQ5D4XBAkUKQIjR54mMPBX1wQghMjWmjZtytSpU2/7TGtNr169CA0NZeDAgTc/DwoKYseOHQDs37+fH374gc6dO99xzebNmxMTE8PkyZNvfrZjxw7WrFmTbAx79+7Fy8uLF1988Y6EBzBx4kR8fX0JCgpi586ddOjQgf379/PFF1/cLBjbsmVLBg0axKZNm0hISMjw95CUjPSSmDNnDitXxrJ+P1AVKMCttXpgKi5sxtTZqw6PuPj57YYNrzB//g7efXe0awMRQmQLic/0khozZgz+/v7MmDGDatWq3XzWN2rUKLp06cKCBQuoWrUqBQsWZObMmQQF3TlZTinFjz/+SP/+/RkzZgw+Pj6ULl2acePGAXD9+nVCQkJuth84cCAvvPAC/fr1Y9KkSXfMQUhISKBbt25Ur14dMEn54sWLuLm53dxyMfF3XFxclr8bkKR3hzVr1rBoUQxX4oCymKQXyq0xcbDt92mgKjR6KJmLONHKlSeZMmUiAwa8TMGCBV0bjBDCjpO2a0pi/fr1qR7XKWwj9fPPP6fr+sWKFWPOnDnJHks6Ehs0aBBTpkyhZMmSFC5cmNOnT7Nw4UJ69eoFQL9+/Rg6dChFixYlICCAwMBAunTpQnBwMMOGDSMsLIylS5eyY8cO6tati5tb1m9OqpS+gJwiPDxcb9q0ydFXZfxfm+m/JJlDscAooClwH3zZAp6qB56eDg4hna5fB6iMn99uMrOWRwjhGHv27CHUGaVXcrDFixfj5eVF8+bN09V+6tSpFCxYkHbt2qXYJrnvXSm1WWsdnlx7GeklK8+tZ3qJ4jE5xQsz+jsN7IVnRsA9m8BVs3jNJJq9aD2f6OjW+Pomt8BQCCFcr02bNhlqn3RxvSPIRJYkVq1aRa9eh5k8BEjcqGAP8A5w1va+CCbpFYVmvaFo0WQu5ERaQ7NmT9KvXz/XBiKEENmcjPSSOHToEAsXnsHDFzORBUxBWY2ZzFLE9rMH8IOijaFYMdfEmkgpaNPmCoULyy4tQgiRGhnpJdGrVy9On36cX9djqiwA5LP9vmz7XcT2+ywcPQlnzjgzwuQNHgw9e25zdRhCCJGtSdJLVh7y+9z29lYxWbiV9E7D+tfh9dedGlyK4uN/Y/bs0Zw7d87VoQghRLYkSS+Jffv20bPnKoYNAFbYPrQvJgtm5OcJnAaPNvDkk86PMzn790PnzkOluroQQqRAnuklcfnyZX777RjVq4O6brfSpi5mxAcmCRYGTkNMW6hR1xWR3ik0FP78U1G/fltXhyKEENmSjPSSqFOnDsePv82iRZC/o92BhkB1u/eJMzijYNUGp4aYqoYNNW5uH7g6DCGEyJYk6SXLH+D253rxmNubiRsOFAGigNXwUBOIjXVqgKn6+efpNGpUJ9mtiIQQIjeTpJfE1atXeeKJbxg4EM5+YXdgK/AxcMX2PnEySxC8ONapIabJxyeehIRITp065epQhBAiW5Gkl4Sbmxvr1h3k/Hlwi+fWyM6+xBCYZ3oAMVC8Hnh5OTXMVLVoAX/+eYlSpQLTbiyEEMlYvHgxlSpVonz58rz33nsZble6dOmbG1uHhye7I5hLSNJLws/Pj0OHZjNtGrR4m1vfUNKk5wcEAKdg21bIToMqpUCp61y7NpHjx4+7OhwhRA4THx9P3759+fXXX9m9ezczZ85k9+7dGW63cuVKtm3bhuP3R848SXrJMtM083nbfZTX9vuS3WdFgFMw8wXIQF1Fp9Aaatd+i759n3d1KEIIJ/r7779p0aIFFStW5O233+aFF15g48aNGbrGhg0bKF++PGXLlsXLy4vOnTszf/78TLfLTmTJQjL69PkAf3/4dTHQBJPcfGw/SZPeYaj4HDz4oAsCTYVS8PbbcRQt2tDVoQghnCQ6OpqOHTsyd+5cypYtS+XKlalduzZ16tS52aZx48ZcuXLljnM//PDDm4VeIyMjKVGixM1jISEhyZYsSq2dUopWrVqhlOK5557j2WefddifMysk6SXj8OF/CQkBb29uPdMDaAHYb29ZxBy/mg9KlXJmhOnTsSPAMmCIiyMRIpfp3x+2bXPsNWvWBFux1pQsW7aMe+65hypVqgAQGxvLyy+/fFublKqcO9off/xB8eLFOXPmDC1btqRy5co0adLEKX2nRm5vJmPJkvl89RUMmgzYV1AIB8rYvbdNZjm9HdKo2+gyp0+vYOTIfly9etXVoQghLLZt2zbuueceAE6ePEmePHm49957b2vTuHFjatasecfPsmXLbrYpXrz4bfMBTpw4cUfV87TaJf4uXLgwHTp0YMOG7LGgWUZ6yTLP9Ir4J/n4OnAeSBzNFwTcIH4jdO8Be+58zutyhw/DW29N4J57WvNgdrsHK8TdKo0RmVW8vLyIjIwEYMiQIcQms4A4PSO9OnXqcODAAY4cOULx4sWZNWsW3333XbrbXbt2jYSEBAICArh27RpLly5l+PDhWf8DOoCM9JLx9tujePFFd0b2BP60O7AF+ApIXPPtgUl8/vDWOOfGmF7168ORIz48+GBjV4cihLBY165dWb16NZUqVaJGjRo0aNCA/v37Z/g6Hh4efPbZZ7Ru3ZrQ0FA6dep085YpQEREBCdPnkyx3enTp2nUqBE1atSgbt26PPDAAxkuIGsVGekl4+LFi5w750HxYvFssd/UxH7ZQuI6vSLAMchX1qkhpptSUKpUNDAFrV9GKeXqkIQQFgkJCWHz5s0OuVZERAQRERHJHvvll19SbVe2bFm2b9/ukDgcTUZ6yRg7dizffRfMrJmY53iJkq7VA5P0LsOPc+HiRWdFmHFDhrxD586PuToMIYRwKUl6KcqDnyfksd9pJaWkB3wxFHbtclJomZA37yUKFLhCQkJC2o2FEOIuJbc3kzF16lS+/fY4TZpA3FdAD9uBPIDi9qRnu81Z82GwTZrKloYMAYhF/p0jhMjN5G/AZMTHxxMT40axYhBYwu6AO/Aot5cYygv4wMXL4J90tme2s4IDB37h2rVrrg5EiLuS1jrtRsJhMvN9S9JLRq9evVi9+j569YJ7+yU5WAUoZPdeAUXg1C5Yt855MWbG7t1QseIDUlldCAv4+Phw/vx5SXxOorXm/Pnz+Pj4pN3YjtzeTJEZtt2xVu+c7aey3WdFIHoDfPopNGjgpPAyITQUPv/ciwcfbO7qUIS464SEhHDixAnOnj3r6lByDR8fH0JCQjJ0jiS9ZKxcuZI33ljOq6/Cd88CDwClbQe3AWuBN7g1TrZNZsnEchinUgqeey4WWAFUcnU4QtxVPD09KVOmTNoNhUvJ7c1keHp64uvrTf78EFYPs9F0okDMfpz2u3rZJrPsPuy0ELNkzZqP+OKLL9JuKIQQdxmnJj2lVBul1D6l1EGl1OAU2nRSSu1WSv2tlLpz3xsnaNSoEcuWdee+++CVMUCw3cHEEkOX7T6zJb2p3zgnvqyaMeMgb789nLi4OFeHIoQQTuW0pKeUcgcmAG2BMKCLUiosSZsKmJIA92qtqwD9nRXfnWz7b+ZJ8nHiM77rdp95A77w+1JIpmJHtjNqFBw40BZPT09XhyKEEE7lzJFeXeCg1vqw1joWmAU8lKTNM8AErfV/AFrrM06M76Zjx45Rt+4Efv0V+rQHFtod9LP9vp7kpBDIVywnLFuAggXB13cekAMytBBCOJAzk15x4Ljd+xO2z+xVBCoqpf5USv2llHLJDqU+Pj4EBeXFywvatQPsJwflBboDFZKcVAwu/gsrtsP1pAkxG9q//xoNG9ZMtjCkEELcrbLbRBYPTDppCnQBvlRK5UvaSCn1rFJqk1JqkxXTg4sUKcKvv77F/ffD6HfAx37/TQ/MTM6kIzpbEmwZDjVqwqVLZGtFi0Jc3L9cyu6BCiGEAzkz6UVyqxIdmPFTZJI2J4AFWus4rfURYD93jqnQWk/WWodrrcMLFSqU9LCD3HqYV9g3yaE9QNKZmiHAC0AROHgAypWDTz6BZMpZZQsBAbBxYxStWuV3dShCCOE0zkx6G4EKSqkySikvoDOwIEmbnzCjPJRSBTG3O12yEKBu3f589BGMGAEnhgL2myysBJIrAlwAeA68+0DFMHjpJbMgfPZsyK6bNMTHf35b5WMhhLibOS3paa1vAP2AJZix0hyt9d9KqZFKqfa2ZkuA80qp3ZjU8qrW+ryzYrRXvnxpgoLMDivl2mLW5iXy5c6JLHZiisC1evDII5AnD3TuDI0bg4PKXDlUhw7TaNcuQrZOEkLkCk7dkUVr/QvwS5LPhtu91sBA249Lfffdp0B5AO6LhgNb7Q76YbYiS8WOTXD6MPy9BX76CYYOhTp1oEcPePddCA5O/Xxn6d07nqioBmitpcCsEOKul90msmQjt57pFfIF4u0O+ZHqSA+ARnD9aYjyhF69YP9+ePllmDEDKlaEDz6AmJg0ruEEERHwyCObcHOT/xSEEHc/+ZsuBZ0796VnT1MY9r3WwF67g4lJL7U7gh5wJQGe+QkWLYLAQJPodu2CJk3gtdegalX4+WfXP++7fHkrkyYNkpmcQoi7niS9FFSqFEaFCoqQEOjwLFDQ7mBdzEzNdFj8rVnrl1hVvWJFWLgQfvkF3N2hfXto2xb27HHwHyAD9u+H559/n59//tl1QQghhBOonD6BITw8XG/atMmiqwcAV1l1FJpNy+QlYqDIeTg2AbyTPEGNi4MJE+Ctt+DqVejXD958E/I7eRWB1rBrly9Vq55BqaT7rgkhRM6ilNqstQ5P7piM9FJlEkABLyDa7uNLwJ/AxXRcwhtOF4PPN915G9PT05QjOnAAnn7arOurWBE+/xxu3HBE/OmjFFSrFoVSc53XqRBCuIAkvRS8/vrr1K59AYAurTErCBNdAX4DMrAz6Bvj4b6mEB9/57FChUyi27IFqlSBPn3M874ff3Tu877hw9/k7bffdl6HQgjhZJL0UlCpUiUaNTJ1hF57BdzvsTuY0qbTqbh6A05cgwsXUm5TsyasXAnz55vR18MPQ6NGsHZtRqPPnEOHjnPkyHbndCaEEC4gSS8FTz75JOPHm+riTz0FwfZ3hzOR9KgCpx+GOJ/UmyllJrfs3AmTJ8ORI3DvvdChg7kNaqUZM+Drryta24kQQriQJL1UmWd6sbGQ377eqjfmm8tI0lNwPQ4G/2xGc2nx8IBnnjGJ7p13YPlyk/xOnsxAnxlklupN59q1y2k1FUKIHEmSXgqmTZtGoUIrOXcOPv4Ydg0GEjePVpityKIyft3/jYIOj0BUOs/194fXX4f16+HaNeja1dpJLrNmRVKoUGFOnDhhXSdCCOEikvRSUKZMGTp2LI27O7RpAw36JGnQG8hEtT/dDGoNBt+klRvSEBoKkybB77+bJQ5WqVMHnnkmRLYkE0LclWSdXqr6AhMBGLwMxvzpuCv/0QMaljDP8DKiVy/45hv49Vdo3dpx8dzOC1P1qWBaDYUQItuRdXqZZp7pxceD52Vuv525G1iT+Ss/1huefTbj5336qVnW8MQTEJm0GqHDxLJr1wccsHrmjBBCOJkkvRSsX7+efPk+YcUK+PdfeOcR4G+7BoeAvzJ//circPwyJCSk3daenx/MnWueCVr1fC86Gu6990NGjRrl+IsLIYQLSdJLQXBwME89VZdixaBIERgwCiht18APM/LL7N3h++Fyq8QZkxlTubJZzL56tTXP93x84PvvE3j//a6Ov7gQQrhQppKeUqq57XdRpZS7Y0PKHkqVKsX48b2pXNlsF9ajJ7c/4vLFFJbNQnmgdSfg21Vw+nTGz33iCfN8b9QoWLIk8zGkpEULKFToR8dfWAghXCizI702SqkQ4HPgYwfGk80Uunn7Mf4SYF/DPTML1JO6Dk+2gvffz9zpn3wClSrBq69mIYZU/P77DIYMsejiQgjhAplNevmAQcBr3L4V810lJORx+vc3r/s/CyywO5iY9DKxVs/+GgkPQesnM3m6Hzz/vNm9JbF0kSNt3HiVr7+ezIXU9k4TQogcJLNJbyQwX2u9D3OT767Up08P7r/fvH7rLcjf1u5gOWAYUDyLnVSFL/dn/vROncxzwZkzsxhHMp5/Hk6cCKdAgQKOv7gQQrhAltbpKaX6a63HOS6cjLN2nV4cZs8x8x3V+Bx2ZOL5W1rUOWh/BmZMgYCAjJ/fqhUcOgQHD2Z83V/a3IAjaF1CFqwLIXIEK9fpDczi+dma1h5ERQUCcP48eJ8EEksD3QAWAvsc0E80/LoIduzI3Pldu8Lhw7BhQ9ZjSerw4QRq1arP4sWLHX9xIYRwsqwmvbv6n/59+/alTJkrgKltt/Ft4KrtoDuwGXDEFpUhkDAASlXL3OkdOoC3N3z3nQNiSSIkBAoUuCijPCHEXSGrSS9n72GWhgcffJBXXy0BmC2/Hn4Hs1QBTLr3I2uzN+3c8IAP15oRZUYFBkJEBMyenXyR2qzw8oJly6Jo0yaPYy8shBAukGbSU0pdUUpdTubnClDMCTG6TNu2bXn5ZVM9tkQJqN8Usy1lIgcmPYAJb0F43cztstK1q1nvl56yRZkRF/c1x44ds+biQgjhJOkZ6dUCIrTWee1/gLZAZWvDcy2tNZcuBRIbC3FxcOkA8J9dAwcnvRsVoGLbjG9NBvDAA2YSjBWzOAHatJlOx46PWHNxIYRwkvQkvY+BS8l8fom7emE6LF++nHz5prJ+vRl9vdsDsF8P5+g7fuVhYwmIycRNY19f82zv++8hJgu7xKRk4MB43nzzPnJ6VQ4hRO6WnqRXRGu9M+mHts9KOzyibCQsLIwPP3yIUqVMUpk0C6hu16Aj0MOxff53HV7+FLZsyfi5XbrApUum7JCjPfAAPPDANpnQIoTI0dKT9PKlciyDpVBzlmLFivHyy49RsqR53z4CCLS403j46i1TQiij7r8fChWy7hbn2bMrmDDhbW5YWbpdCCEslJ6kt0kp9UzSD5VST2Mm7d/V/vvPh8RduE4eAI7YHTwIzMSxG7F5QsJT0P7lTJzqCR07woIFcOWKA2Oy+fNP6NdvOH/+6cBqukII4UTpSXr9gR5KqVVKqbG2n9+BXsBLlkaXDZQp8xRvvmlej34X3O1vHV7FLE534GQWAArDhEzc3gRzizM6GubPd2xIYJZF7NpVgvvuu8/xFxdCCCdIM+lprU9rrRsCI4Cjtp8RWusGWutT1obneh9/PILHHjOv330XSj9nd9ARlRZSsHw51AjP+IitYUMoWdKaW5xeXlClynGyVDJeCCFcKN2L07XWK7XWn9p+VlgZVHbSo8cLNGpkJm9Urgwly9kdtDDp4QWRF+FEBnd8cXODzp1h6VI4d87xYcXGwksv9Wb69OmOv7gQQlhMKqen4dKlaxw+bHaBjoyE6G3cKhybOI3HiqQXAle7Q6FSGT+1a1ezxGLePIdHhZcXrF+/lwMH/nb8xYUQwmKS9NLw6quv0qDBNQDWr4d144DE8nJ+mNmcFs3ij4mHievg+PGMnVe9OoSGWrMXJ8DatQm8/XYVay4uhBAWkqSXhh49evDZZxUAsyTgpelAIdtBX2AAUMO6/t/pDt0zuBZQKTOhZc2ajCfM9HBzA5hKdPRdWz9YCHGXkqSXhgYNGtCxYyhgNnauXhXwcF7/cQ2h3mMZP69LF/N79mzHxpPovfdWUr58GeLi4qzpQAghLCBJLw3Xr19n1y53rl0ze2LuWQFE2jX4GfjNwgDCYEUmbp+WLw/h4dYtVA8PhyeeKEdUVJQ1HQghhAUk6aVh7dq1VKs2jy1bzG3DT9/g9v03z+GYmnqpWH8Q+gzO+GzMLl3Mdmb79zs+phYt4L33IsmbNxOl3oUQwkUk6aWhRo0azJnTncqVTdL7bR3Q1K6BL9bM3rR3GT4fAxktXv7YYyZmq0Z7Wh9l06YvuJC4ZY0QQmRzkvTSUKhQITp2bEUh2+SVulUBb7sGDi4vlKwi4D4AmrXP2GnFi0OTJibpWVEcYf9+qFOnj6zZE0LkGJL00mH79iscse25ue4P8LW/vekHRGF5Dfn4QJi8OePJq0sX2LcPtm1zfEyVKsHs2V507/6w4y8uhBAWkKSXDk2bvsJHH5nX330HsUvsDgZh6sc7YRLjR2OgQwbzy6OPgoeHdbc4O3WKJV++pdZcXAghHEySXjrMnPk5ffqY1++9B3VH2x28B3ga8LI+jqsa/kswu62kV1AQtG4Ns2ZlriJ7evz001gmT55szcWFEMKBJOmlQ5s2nQgLM+sGChSA4oVdFEhDoIMZuWVEly5mkfratZZExaxZe5k48WOpqi6EyPYk6aXDoUPH+OOPPAAcPQqRvwKXbQfPAhOBw86JZfUxWLoZrl1L/zkPPWQqv1t1i3PiRNi8+SGpqi6EyPacmvSUUm2UUvuUUgeVUoNTafeIUkorpcKdGV9K3n//fR55xEzRPHoU1n2DWZ8H5hs8w60kaLWz0KYu/O9/6T8lTx548EGYOzdjt0bTq0ABcHf/HxDv+IsLIYQDOS3pKaXcgQlAWyAM6KKUCkumXQCmOO16Z8WWlv79+/PTT2aD5XvvhY9WAWVsBxPLCzlrY5KC4NkWGjXP2GldusDZs6ZOnxVWr46kWrXynDx50poOhBDCAZw50qsLHNRaH9ZaxwKzgIeSafc2MAbINrsZh4aG0qCBKaTn6QklCnGrsoKP7bXVa/USKYitAysyuB68bVuzd6hVtziDgyFv3iucPXvWmg6EEMIBnJn0igP2e/6fsH12k1KqFlBCa73IiXGl6fz58/zySzT//WfeL5sOHLAdVDhngXoSY783z9LSy9sbHn4YfvgBrNgus2JF+PPPq9SoUdLxFxdCCAfJNhNZlFJuwEfAy+lo+6xSapNSapMzRhZbt27lgQd+ZedO837OFG4lPYDyQAHLw7jNsTUw+HW4noFk26ULXLkCv/xiVVQxREVN48yZM1Z1IIQQWeLMpBcJlLB7H8Lt9QoCgKrAKqXUUaA+sCC5ySxa68la63CtdXihxP3BLFSnTh3WrRvAPfeY9zv2AxF2DToA91oexu2aQNNx4OeXZsubmjWDIkWsKy574wZUrPgaQ4cOtaYDIYTIImcmvY1ABaVUGaWUF9AZWJB4UGt9SWtdUGtdWmtdGvgLaK+13uTEGJMVGBhI/fp1CbAVFCiWH9xcPTvfD345BicyMGvUwwM6d4aFC+HiRceH5OEBw4fH8eST9R1/cSGEcACnJT2t9Q2gH7AE2APM0Vr/rZQaqZTK4FbKzrdw4SG2bjWvf/oRfOwXei8HPnN+TPGXoVET+P779J/zxBMQGwvz5lkT0zPPQJMmW625uBBCZJFTn+lprX/RWlfUWpfTWr9r+2y41npBMm2bZodRXqJu3cbw1Vfm9bJlEGcfmQYuYPmm03fwh5NXIDo2/afUrm02is7IOr+MOnnyf3z11STrOhBCiEzKNhNZsrvff/+JYcPM688+g4Yf2h30AxKAGCcH5QZxT4Jb1fSfopQZ7f3+O/zzjzVhzZ59maeffp7Dh520TY0QQqSTJL10ql69CUWKmAd5bm5QyN/uYOJkEicvW0g0YYMpH5ReXbua31ZNaOnRAw4cqEPZsmWt6UAIITJJkl46rV+/mblzTab7+2/YOw24aDvo4qT351dQO9wsR0iPsmXNzjL/+581xWXz5YPy5Tdy+7oOIYRwPUl66fT111/Tr59Z1X32LBxcDSQmmfxANZxSXihZNaDec2YBeno9/rhJ3jt2WBPS5cvw7LOP8tNPP1nTgRBCZIIkvXR666232Ly5FgBNm8KHy7i16rAQ8AjgqpJDwbC+AERloF5ep05miYFVE1ry5IE///ybw4dltCeEyD4k6aVT0aJFCQm5tbb+tmd6iVxYTu5aFPR5Bzalc75rUBBERJjnevEWFEdwc4OdO+MZOLCC4y8uhBCZJEkvnY4fP86kSWc4c8Y8B/vuPWCX7aAG3gNWuC4+EmD2xxmbnPLEE3DyJKxaZU1Ibm4AU7iWkeJ/QghhIUl66bRv3z6ef/4P9u410/63reVWTT0FuOOyiSwAeEHCMxDRN/2ntGsHefNau2Zv2LBfCAurRLwVw0khhMggSXrpdO+99xIZ+TYNG5r367cDTe0auKDSwh3yw6TN6W/u6wuPPmp2dMnIxtUZ0bixpmfPisTEOHsRoxBC3EmSXjr5+vpSrFh5PDzM+yC/WyX1TANcn/SAH+dCrTrpr5D+xBNmqcPPP1sTT6tW8OabR/Hz87WmAyGEyABJehkwceK6m8+/pkwGr6V2B7PDSA/QnnAJs6wiPe67D0JCrL3FqfUR1qwZJwVmhRAuJ0kvA954Yypz55rXR46Am33ZuMqYwkiuVgmud4aC6Vw+4eZmdmhZvBjOnUu7fWYcPgxNmgxkypQp1nQghBDpJEkvAw4cWMf48eb1mDEQbl82riZwnwuCSsapqzB9PaS3luvjj5vbobNmWRNPuXKwaJEbL73UyZoOhBAinSTpZUBQUEU8PG49ybttrZ4GojEbT7taHPRuBW+/nb7m1atDzZowdap1IUVEJODn9611HQghRDpI0suARYuW8MknZkLGhg2weSzwn+3gLsxaveywAYkn3GgOjR5K/yk9e8LmzbB9u3VhLVgwnldffdm6DoQQIg2S9DLg559/ZvRoM/U+NhaizwOJtewqA0HArzi/xFByasOSDKwJ79oVvLzgm2+sC2nbtgssWDBLFqsLIVxGaSu22Xei8PBwvSm9e29lUVRUFN7erXBz+wOAT9bDS4vtGhwDvgHqAW2dElKqvK/AAD94963E3VFS99hjsHy52aXFy4LNs6Ojwcvrftzcljn+4kIIYaOU2qy1Dk/umIz0MsDX1xc3t1vTIgv5JWlQCqgDrAeOOzGwFMQcg/dHwbZt6WvfowecP2/dmj0fH3BzW0F8/B5upHchoRBCOJAkvQzYt28fQ4YcITISEhJg/EvAliSN7sfc5rzg/PjuEApFXofqNdPXvGVLKF4cvv7aupAiIzXly9fn229lUosQwvkk6WXAiRMn+PDD7Rw7Zm4X3ogCkg5YfIDngRrOj+8O7vCvG/ywJ33FYt3doXt3s2YvMtKakIoVg/vuiyYkpJA1HQghRCok6WVAs2bNiI39+Ob+mz8vBeom09Ads4RhJ3DaaeElT8MLz8GAAelr3r27GcVOn25NOErB1Kmx3H//KWs6EEKIVEjSywA3NzeUKnLzfcGkz/TsxWBmci7AtWv3FJyJh2sq7aYA5cubrcm+/jp9o8PMunLlM1ZZVdNICCFSIEkvA7TWDB78E/Pnm/ej3wXPH1Jo7IOZwRmJmdjiSm0gNgO7xfTsCQcPwh9/WBfSoEFbiYhow6VLl6zrRAghkpCklwFKKaZPX856WxLz8ACf1Kb2VwUqYorLXrQ8vFTN3gUr16evSvojj0BAgLUTWgYOhBUrWpA3b17rOhFCiCQk6WVQZOQORo0yr4cOhWp9UmmsgAjM873frI8tNTEHoHl9WLQo7bb+/tC5M8yZY8oOWaF8eahffwVKXbSmAyGESIYkvQxSqiD2lfTuWKuXVD6gJRBqXUzpUgoCO0CDe9PXvGdPU1g2saqEFa5cieL11x9lxYoV1nUihBB2JOll0PTp3zF0qDdgdi9ZM5S01+TVw/Vlh9zhUg1YdjJ9zevVg9BQa29xenvD9Om/s3athQ8PhRDCjiS9DNq6dStLbcVj8+SBoGDM7cu0xANrgG2WhZYuI7/hZnmk1Chldmj580/Yt8+aWLy8YM+eeN54o4I1HQghRBKS9DLo448/ZtOmOoAZDT3/EWYHlrS4YSowLAWirIsvLXvXwPsfQVxc2m27dTOTdb74wrp48uQB+Eg2oRZCOIUkvUxJZf/NlCjMEobrwO8WhJRe90Otd8HTM+2mwcFmJuc335jne1b59ttNFCtWhFOnZMG6EMJakvQyaOvWrTzxxA6OHTPVxt/oCPyVzpOLArWBDcBZy0JMnS8sOgR7zqZvtNe3L1y8CN99Z11I9epB166FyOkVP4QQ2Z8kvQy6evUq69ad49w5c+uvUiiQJwMXaA54AovTamgdHQX31oEJE9Ju26iRqaw+YYJ1O7SULw+TJv1D0aJyi1MIYS1JehnUuHFjDh16n9q1zfsp08nYzEx/oD3QzPGxpZsvXCoIeQun3VQpM9rbtg3WrrUyqAT27BnG6tWrrexECJHLSdLLlOI3XxX2z8TpVYAQhwWTKQkPwoF0JD2Axx+HwMD0jQyz4qmn5vDSS/2s7UQIkatJ0suE3r2/uTntf8CL4D4jExe5gdmM2oX7ck5aD9/OTvu2pb+/Wb4wbx5YOdfkq68SWLr0Qes6EELkepL0MuGff/7jzBnzukoVCCiTiYt4AJeAlZgZnS5waTM80RnSc0fx+efNxJcvv7QunmrVoFChb4BY6zoRQuRqkvQy4ZdffuPdd30Bkwwqd87khVoC0ZhF665QBQo+k76tySpUgNat4fPP0zfrM7OOHfuXhx6qz5YtSUvSCyFE1knSy7Rbz/XSvVYvqWDgHswtzrS2MrOCB5wrDrP/Tl/zvn3h5ElullayQmAg7Nz5N0ePHrWuEyFEriVJLxPmzZvHffedJS4OfvoJljwLnM/kxZphKq27cM/lwe/Byy+n3S4iAkqXtnZCS758cPBgLA8/HGBdJ0KIXEuSXiYkJCQA3ly5AqVKQdXmmGd0mZEXeBho4bDwMuzkMVi1Ie1ae+7u0KcPrFoFu3ZZF4+bG2g9kv3791vXiRAiV5KklwmdOnXi99+7U6AA3HMPdB0EBGbhgqGYEkSa9G1e7WgtIe/TJqmlpVcv8PGBiROtDemjj/6gSpUwjh07Zm1HQohcRZJept1aaFfYn6wnq2vADGB3Fq+TGW6w6igs3A7nzqXeNCjIFJidPh0uXbIupE6dYNy4MhQpUsS6ToQQuY4kvUy4du0aTZpMYto0iImB3o2ArJaE8wWuYiqs38hyiBkXCw83hmHD0m7arx9cuwZTplgXTokS0LfvQXx8XLiQUQhx15Gklwl+fn54efnj4WEKoT78OFAsixd1A1oBFzEbUjubF8Q1hdrpWBteuzY0bw5jx5qkb6UFC15g7Nix1nYihMg1JOllglKKZct+4vHHzft3xwDlHHDh8raf1Zjbnc5WG2acSV/ToUPh339h2jRrQ1q4cCdTp07kxg1XDH+FEHcbSXqZFoxZa2Bbp5fGzMd0awXEAJZu7pyy1bvh0V6Q1jK55s2hTh14/31TYskqH34IW7eWxsMjs9NjhRDiFqcmPaVUG6XUPqXUQaXU4GSOD1RK7VZK7VBKLVdKlXJmfBkxcuS7NGtmkt4TncHtKwdduDDwGNDEQdfLqBvw47dmWUJqlDKjvUOHYO5c68LJmxc8PFZw48Zqrl69al1HQohcwWlJTynlDkzA1A8PA7oopcKSNNsKhGutqwPzgPedFV9GFSlShLJl8wLw6KOQr6EDL14Z8MZxo8eMCISE/hCSjqTbvj2EhcHo0dbV2gPz3LBatQjeeOMN6zoRQuQKzhzp1QUOaq0Pa61jgVnAQ/YNtNYrtdaJ2y//hcsL8KTsueee46uvGgPQpQuUb+3gDs4Dk4CDDr5uevjCm6vgQhpbo7m5waBBsHMnLFpkXTje3tCt2zVatMi2/zkIIXIIZya94sBxu/cnsN/A8k69gF+TO6CUelYptUkptens2bMODDGjbv0lnN8NSHDgpQMxI70luGTEt/YXKFrM3L5MTZcuZleaUaOsHe0NHQrt2i23rgMhRK6QLSeyKKWeAMKBD5I7rrWerLUO11qHFypUyLnB2fz999+ULfstS5aYhdpLemFKBTmKB2ZSy1nAFQUHSkOBxpAnT+rNPD3h1Vdh3br0lSjKitjYxYwb159//vnH2o6EEHctZya9SKCE3fsQ22e3UUq1AF4H2mutLV4FlnkFCxakYcPK5MsH4eHQpBfg5eBOKgOlMDX3ohx87bQEwKlGsPVK2k179oTChc2zPSudPg1DhnzCzJkzre1ICHHXcmbS2whUUEqVUUp5AZ0xtcNvUkrdA3yBSXjpXDHmGkWKFOF//3uXevXMZI52PQB/B3eigDaYIrMuKi838Ft4773U2/j6woABsGQJbN5sXSwlSsDOnZrXXitpXSdCiLua05Ke1voG0A/zlGoPMEdr/bdSaqRSqr2t2QdAHmCuUmqbUmpBCpfLJswjSa0hj8asr3O0okB3oIEF106HPb/DyHfh1KnU2/XpY5YXpJUgs6p8eVBqCBcvnkJb+RBRCHFXcuozPa31L1rrilrrclrrd22fDddaL7C9bqG1LqK1rmn7aZ/6FV3roYcG0L692Yfy+caYsawVSmP+l3LFzd6GEPAKBBRIvVlgoNmT8/vvYbfFm2bv2nWMMmXK8P3331vbkRDirpMtJ7LkFM2ateD++33x94eXhgNlLOzsX2AccMDCPpLjCWeA0X+YZ2qpGTAAAgJgyBBrQwoNha5d4wkNLWxtR0KIu47K6beIwsPD9aZNm1wYQXVgJ0cvQpnxFnZzA/gcs3zhecDTwr6S4b4cCuyHQwdMYkvJ6NFmecGqVXDffVZH9RzmSxFCiFuUUpu11uHJHZORXhYlJBRDa/CKAy5b2JEHEAH8R9bLGGVCfCUo1gq80pih2r8/hITAK69AgiPXLSbj/PkveeWV7rh2raYQIieRpJcF33//PT4+v3HgAPR6EtRsizssC1TDJL3zFveVVAhsLweb05hT6+sL77wDmzbBbIu/jzNnEvjssxksXy6L1oUQ6SNJLwsqV67Myy83MM/0XoIgR29FlpxWmFGfKyqsAz0+SvuZ3RNPQI0a5janlfX2QkPh2LEEOnfOZ10nQoi7iiS9LKhSpQqjRz9F8eLQpg2UqueETgMwz/QaO6GvZOzfCl9/C5dS2X3G3R0++MCUJ/rsM2vjKVIE4BUOHtwnSxiEEGmSpJdFCQnFiIqCq1fB9yLO2Scz0Pb7HBDthP7sNQKPfuDpl3qzli2hdWtzqzOtjauzat26vwkNDWPWrFnWdiSEyPEk6WVR0aJP8uqrZn3aH4OANBZxO8wVTBWGVU7qL5EnnIyCd1fBtm2pN33/fTMifPdda0OqWxeGD3ejVauy1nYkhMjxJOll0eDBA4mIgIcfhnGfmmriThEA1ATW47xEa+e9YdCoMaQ2cbJ6dejRw9ziPHzYuljc3WHYsBsEBb2K1hZPGRVC5GiS9LJowIDXiYjwISAAXuoHvz4OFTyANBZyO8T9gC9mB1Mnlx9KqA+Fu0D+oNTbjRxpktLQodbHdOHCGtq2rSo7tQghUiRJL4sSEhI4c6bIzfcFfCH/IvCYh/WJyA+zdu8kzl+7lx+OFIf3/4QbN1JuVrw4vPyyWb7w++/WhhQQAFFR+7h69ai1HQkhcixJelk0YsQIgoOPERdn3isF076GebOhSF4nBFAVsymMCwrNArw5BcqUg8g7ikTdMngwlC0LvXrB9evWxeLpCatWJfDUU0us60QIkaNJ0suiBx54gE8/DSfeLulUrgwPNYOl3cDvEHDN4iA6AM0t7iMFNwLhsj9EpTKL1N8fpkwxVdiHDbM2HqUAfmPx4oGsWLHC2s6EEDmOJL0sqlu3Ln37NsfH585jhTXEzwUPq289Ktvvwzj/NmcQXH4MfkhjMk2zZtC7N3z8samybqUbN+Dll8fz/vtvW9uRECLHkaSXRQkJCZw6FcDFi3ceCw6GlStg3mTwcMY3vRtYBhx1Ql9JDF8CPZ6Hf/9Nuc3775tCsD17QrSF6ws9PGDhwgTmz89jXSdCiBxJkl4WXbhwgaJFhzFtWvLHGzSAh6rCpFbAX4CVM+pbAfmBH3H6ovWY/2D6V7A4lcdpAQEweTLs3WtmdVqpTBnw9l5ITMyn/PXXX9Z2JoTIMSTpZVFQUBATJrzC/fen3s57P6glwHELg/ECHsZUe3D2XI6CkPAinK2QerPWrc1I7/33YfNm68N67bX+3H9/M06lVfpdCJErSD09hzgBlEi1hdawdSvMOAPj1lscznJgDdAdU3XdiTzcYFJtqFEo5YX6Fy9ClSoQFGSqMaRVrigrTp2C9esL8tBDe4CC1nUkhMg2pJ6exc6d82TfvtS/SqWgVi34qDW09gN2WhjQfZgZnSUt7CMFN25A757w3PMm0ScnXz744gvYuRNGjbI2nuBgeOihc0BXjh07QoLVRf6EENmaJD0HePHFAUREqLQbYpKf1x/g+wemGroVPIAamP91L+Dc53tuEP8IqC4Ql0p+adcOHn/c7Mv5hxNmnB448BtVq1Zm3Lhx1ncmhMi2JOk5wAsvvMBnn6V/s+P/zYDtf0G41SOxGOArzMQWZw5wCsOWK9BvEaxalXKzzz4zE04efRROnLA2pPLl4bXXYunUKb+1HQkhsjV5pucwD2OyS/qduwZVusCZokAxS4IyM0YXY255NrOoj5RsBeab7ceaNEm+yZ49UK8eVKoEq1ebyuvWyofWG7l2LZg8eWRJgxB3I3mmZ7Hr16+zcaNHsmv1UuMeC57bIN8BK6KyqYepxvA7sMfCfpJTHTweBp9UBsGhofC//5kJLb17p/wc0HEu8uKL9WnUqCHXrlm9VY4QIruRpOcA27dvp27duaxdm7Hz8ueHzRthyzwIsWqfTgU8ABTHDERTKQXkcO5wozp0nAc7DqRchqh9exgxAqZPh/HjrQ/rwQfP0779Zfz85D9/IXIbub3pAJcvX2bVqreoX/9jChfO3DXW7oFmHSG2FbcqozvSZcwWZS0w6/mcKQF8p0CtsrBmdeL+mEmaJJhnewsWwJIlpLnu0TEe4NSpSeTLVwif5PaRE0LkSHJ702J58+alffsHM53wAHyjIe9ZyHPVcXHdJi+mDJEXZgPsGIv6SY4bRN0P15tBVAozVt3cYNo0s1l3p05w5Ij1YUVFLaJRozB69eppfWdCiGxBkp6DbNt2ma1bM3/+PffAP0dh2XAI8MK6UkEJwAzgW5yb+MrBVnd4aBbMXwjJPU4LCICffjKjvvbt4fx5a0Py9YVBg67y0ktXgZx9x0MIkT6S9BzkmWdGMmRI1q7h6wv1QuCtYuD+BXDJIaHdzg1ohNkO7Tsg1oI+UrFsC3T4PxiRwt6b5cvD3Llw4IC5xXnunLXxPPMM1K37MzCQDRs2kNNv9wshUidJz0EmTvycjz5yzGyU+hWgXigEW7VrVlXMCot/cH7iywe6K+yvBjdSWDvYooV5trdvHzRvDmfOWB/Wxo3jqF+/Hl9//bX1nQkhXEaSnoPUqVOHsDDHrDZv2BD+XAXrnoeyeQErRjvVgP/DlCFabMH1U1MO5h+Cx2fD6NEQFXVnk1at4Oef4eBBU4vv9GlrQwoPh08+gS5d/sJlZeiFEJaTpOcgkZGRzJ7tnuxf4JlVOh+0OQZuU4ArjrvuTTWAjkBTC66dDnMWwNDXYeny5I+3aAGLFsHRo9C0aeq1+rJKKejXD/z8phAT05bnnuvBP//8Y12HQgiXkKTnIFu2bKFz5+1s2uTY9QBDXoEPxkCDUIde9pYqmJmdCcBvWJNcU1IZ6AsjT8KJy3D9+p1NmjWDX3+F48dN4ouMtD6svXt/Y86c6axfv9D6zoQQTiVJz0HatWvH77//TuPGe4AIli+Hy5ezft2QEBj4IvzWDVrkAVZjzT6aZ4ANwGRMpSRnKQhb/oWab0NIKViczK3WJk3M2r2TJ6FRI9iyxdqQatSAQ4cS6NhxJLCec1bPphFCOI0kPQdRStGkSROgLJcufUeHDr4MGODvsOv7e0G9/yDPNnC3ompCMNALcAe+weyb6UTngUtFIKUKTffeC8uXm9JFDRuaCuxWTrQsUADgNPv2NaF8+VJMnz7dus6EEE4jO7IA9O8P27Y5Ipybrly5gqeXOz7ep4iNPU5cHPhnMQdqIDYWohJg91mIjcGUEXKkBMxWZdGYnWHyOfj66VAiEHyuQ6HC4OV5+7G4ONi9B/77D4oUgYoVwN3dulgSEuDQIShRshA+3pXQuJO+IlJCiEypWROyWAJMdmRxgYCAAHy8/YCyHD1aiM1bIC6L9fMU4O0F+XyghDsQieOXG7gBRTDP+SyveJC842fhwEE4cfLOY56eUL0alC5tZnRu2ZL8s0BHcXODChXAx/sssJF9e3dw9Ngx6zoUQljK0eOEnMniwqJB58+z988/8WzvDgxgxYoDNGkCHln49v3Ow6rP4Eo9+Gg9ZoRm1T9hfgMCgLoW9pHUOfAuBG81h6ZeEFzYJDowyb80cOA3iOgKUbth4kTo1i35fT0d5caNGN57OobSpcvw1ls/Y74UIUROIrc3nWzPnu1UqXIP773nxWuvOWYfsF93QoeWENMEMxvTkeKBmcBBoCjQDlOxwVk0+H4NJQJg7647k1pkJDz2GPz5p3nuN3481K5tcUgalCrJxo2vMmPGfkaOHEm+fPms7VQIkW5yezMbqVy5Ot9//z29e+8CerN9uwerV2ftmjWCoGlNeLE15PHCsdtIugOPA49iljN8CSzCPPNzBgVRj8DBxvDab/DfVVi48NYkluLFTZHaKVPM1mV16kCvXnDqlIUhKYB/WL78BebP/xJv750AsoWZEDmAjPRcrHPnB1mxYhnHjt3A1zeLD/2AyMvQsgfsOYoZlTnynzXRwApgG9AHyO/Aa6dTwG64MgeWrIBWSSrBX7oE77xjRns+PvDGG/DSS+DtbV08UVFmz1StH6R583959NHu9O3b17oOhRBpkpFeNvb117P59dc/8PU9gNY9GTBAsXlz5q9XPC88HAqty0LdErYPHVUg3AdTnuglbiW8WcBKB/aRhiuVgC7QeQO8/TtM/BJmzjQjv8BA+OAD2LUL7rsPBg2CsmVNIrRq/05f22SfK1d+pmDBTeTJMx3YS3R0NHv2OLtUvRAiLTLSy0aOHTtGnTq1GDWqOk8/vYXY2Mtcvw6ZeVxknjvBvHXwWBNIiABqOTjgOOB7YC9mStQ9QAOggIP7SYX7dChdEP5YDsF5zFZlRYuaY7/9Bh9+CEuXgpeXefb3wgvmFqi1FDNnhtK1627WrVtF/fr3Wd2hEMJOaiM9SXrZzLVr13B3d8fHJ55ZswbSs+cUNm5MoEomJ6icPw8ffQT3PwqzTsC0ZRD7D2bCi2daZ6fTWWAdsB0z8aUrUBHzbNHqRW0JQDS4+0PTYFj9Arw5Al4ffKvJ3r0wYQJMnQpXr0K9eua534MPQnCwNWGdOWNKJPXpkw83tyf49FN//v77IpMmTUJZOcVUCCFJL6favXs306ZNY/ToTri5TeHzz6fxzz9RvPtu5qfmP9cPpn4F9cbD2lMQfwbwAxyxecwVYDNQD7PGbz2wG5NgQ7F+hn80sAV8K0PHZtDYB2Z+CJ9+AmFhZlu4adNMAty3z5xSrx489JApWhsWZt2Sh+HDYf36AJYs6Q9E8P77vxMcXJQnn3zSmg6FyMUk6d0l+vbtw969G1i+vBmwkM8+20fp0tCuXfqvobXZYaR8eTh7De67H/45DjG9bfXtzmJ2YnHEvtnbgD9t1wQoCJTFPBd0hgPgthxavAnt7gG9B1bOhy+/NLdB5883dfs2bjTNy5UzhWvr1TM/oaFmcbqjJN5yBqhf34OwsJJ8/fVIoBWvvDKG5s2bExHhrC9HiLtXtkl6Sqk2wHjMRPgpWuv3khz3BqYDtTHbMT6mtT6a2jVzU9IDSEhIwM3NDa01lSqVpUmTIkyZEgxspG/fk7RqZUYu6bV9u7kVV6cx/PEPPNEEvErC9f+Da3HAfqAQWZupeQbYhylaC2YJBMBszO3PwrY+CmGeB1q1ZcJW8NgEHT+GhqVg07ewcy389AP88gvMng2bN9/aKDwgAOrWNQmwWjWoWNHszhLgoBFrdLSZZRoVBeXKudOvXxhDhz5FTEwNmjQZwhtvDOfBBx/kxo0bHDp0iDJlyuDl5dgqHkLcjVJLek7bkUUp5Q5MAFpi9vHfqJRaoLXebdesF/Cf1rq8UqozMAZ4zFkx5gRutqGHUoo9ew5y5coVIB9RUVH8+mtlSpVqCpQkKuovwsJWMnJkAt26aWJjzaimQQOzti1RjRq3Xj9QAb772syCrNcA1h+BJpWgVifIVxv2noYTnwN1MGWB4oGTQBDmFmlKCtt+7GnMaPIY5hZoojCgk+31Ysxt17xJfjL7LPIeuHEPzNxtfjgGKGj6I9QMhtPxkCcIvlwAp/fBt1/D7t2wahXE29WVLVr0VgIsUQKKFTM/RYua34UKpW+E6ONjfvv6QmRkPPHxO4FXOH8e8uaFhITnga84dCiIypW/Zvr0oXTr1o9jx2IZMWIE/fv3p3r16ly+fJl9+/ZRpUoV/PxS+x9CCOG0kZ5SqgHwlta6te39EACt9Wi7NktsbdYppTyAU0AhnUqQuW2kl5bEkeC5c+d46aWX6N69Ky1blmfv3lWEhvZm+vTWdOsWyN69R2jUaAvTpvnzwANXOXIkgVGjzLq2qlVv3f4rV878BZ8nj1kE/vY70Pk5KF4L/toAE/pBw77gWwuOHYQjU4A2EF8auAisxSTJIsBVzIivAuYWahTmf+EgzJKHxCKx1Wy/x3LnIvgGmH82xQBfY54d+mKWU3hhkmYZ23l/YxKkO+afd+6Y5JsXs2fpf5hFO+623/9iJsZUAeLB/SfwdzN/Pr9rsPID8HCHKrXh3HHYt8NsgJ2Uu7uZcevmBvnzQ6lS5vfly1CwoNlOzd8fzp0zm2aXKGHWEl64YNoVLmxmm167Zr73hARYswbuucfclt6zx5sePW7wxRcVady4MmvWxPPQQwtYvnwIjRrVY/nyw/TqNZqFC2dxzz1NWL9+E+PHj2fMmDGULFmSXbt2sWjRIp599lny58/P4cOH2bx5Mw888AB+fn78+++/HDlyhPDwcLy8vLhw4QLnzp2jXLlyuLu7c+XKFa5evUpwcDBKKaKjo4mNjSUgIAClFDdu3CAhIeHmqDTx/74ygUc4S2ojPbTWTvnB7Okxxe59N+CzJG12ASF27w8BBVO7bu3atbVIW0xMjN66das+e/as1lrro0eP6j59+ujt27drrW/o9et/1cHBhfSqVZO11iv1smVva0D//vvLWuuP9S+/9NCAXreum9Z6gJ4/v60G9Cef1NX//NNBz5nTUAO6VatCeu3ae/WUb6ppQAfm99Af/6+C7j8qRGPGd7r/Z8V0x9cKapR5//jYgrr583m1cjPv272XT9fu4a/dPMz7+gP9dalmnjeP3/uevy7dzuvm+b7BSrv7cvP6dd7304Ubedx8f9uPG7r2OD+dv5Z7sseVO7rWRF+dN8wt+eOe6Npf+Gm/snbH1e3n39PJV3sHqGSPA1qpZOJy4k9q/SultIdH8n92QLu5uWkPjxS+W9Du7u7awyP57xbQHh4e2t095eNeXl5pHndzSzk+b2/vVI/7+PhopVSKx319fVM97ufnl+p3m9Zxf39/l56f1nFfX99U/9vw8fHJ9HE3Nzft7e2d6nE/Pz89duzYLP99B2zSKeQMZ470HgXaaK2ftr3vBtTTWveza7PL1uaE7f0hW5tzSa71LPCs7W0lzPghqwoCUi00efLdpEy+m5TJd5My+W5S5ojvppTWulByB5xZZSESKGH3PsT2WXJtTthubwZiJrTcRms9GVPj22GUUpt0SsPhXE6+m5TJd5My+W5SJt9Nyqz+bpy5DdlGoIJSqoxSygvoDCxI0mYB8JTt9aPACu2soagQQoi7ntNGelrrG0qpfsASzNSBr7XWfyulRmLuvy4AvgJmKKUOAhcwiVEIIYRwCKcWkdVa/wL8kuSz4Xavo4GOzozJjkNvl95l5LtJmXw3KZPvJmXy3aTM0u8mx+/IIoQQQqSXlBYSQgiRa+T6pKeUaqOU2qeUOqiUGpz2GbmHUuprpdQZ21ISYaOUKqGUWqmU2q2U+lsp9ZKrY8oulFI+SqkNSqnttu9mhKtjym6UUu5Kqa1KqYWujiU7UUodVUrtVEptU0pZtuNIrr69adsabT92W6MBXfTtW6PlWkqpJph9VKZrrau6Op7sQilVFCiqtd6ilArA1Jb4P/nvBpTZdsVfa31VKeUJ/AG8pLX+y8WhZRtKqYFAOJBXa52B7eLvbkqpo0B40nXZjpbbR3p1gYNa68Na61hMHfCHXBxTtqG1Xo2ZRSvsaK3/1Vpvsb2+AuwBiqd+Vu5g2xDjqu2tp+0n9/7LOgmlVAjwADDF1bHkVrk96RUHjtu9P4H85SUyQClVGlMzfr2LQ8k2bLfvtmHqa/ymtZbv5pZxwGuYXV7F7TSwVCm12bbrliVye9ITItOUUnmA74H+WuvLro4nu9Bax2uta2J2XaqrlJJb44BSqh1wRmu92dWxZFONtNa1gLZAX9vjFYfL7UkvPVujCXEH2/Oq74FvtdY/uDqe7EhrfRFYCbRxcSjZxb1Ae9uzq1lAc6XU/1wbUvahtY60/T4D/Ih5/ORwuT3ppWdrNCFuY5us8RWwR2v9kavjyU6UUoWUUvlsr30xk8T2ujSobEJrPURrHaK1Lo35u2aF1voJF4eVLSil/G2TwlBK+QOtMFV3HC5XJz2t9Q0gcWu0PcAcrfXfro0q+1BKzQTWAZWUUieUUr1cHVM2cS+mNFZz2/TqbUqpCFcHlU0UBVYqpXZg/lH5m9ZapuaLtBQB/lBKbQc2AIu01out6ChXL1kQQgiRu+TqkZ4QQojcRZKeEEKIXEOSnhBCiFxDkp4QQohcQ5KeEEKIXEOSnhBCiFxDkp4QQohcQ5KeEDmIUspXKfW7rSwWSql42+L4XUqpuUopP9vnIUqp+UqpA0qpQ0qp8bZdhxKvc9X220sptVop5eGaP5EQziVJT4icpSfwg9Y63vY+Smtd01bvMBbobdsm7QfgJ611BaAikAd4N+nFbCW1lgOPOSV6IVxMkp4Q2YRSapVSqrLtdVAKFesfB+ancIk1QHmgORCttf4GTNUDYADQM3EkmMRPtusKcdeTpCdE9lEe2G97XR3YaX/QdnuyrNb6aNITbbcn29rOqYKp5n6TrfTRP7Y+ktoF1Mli7ELkCJL0hMgGlFKlgEitdWJx0erAjiTNCgIXk3zmayvYugmT1L7KaN+2kWBs4i73QtzN5OG1ENlDDW5PcrWB2UnaRAE+ST+zFWy9SSm1G3g0yWd5gZLAwRT69waiMxayEDmPjPSEyB5qYktoSqkKwEMkub2ptf4PcFdKJU18SS0H/JRST9qu5w6MBaZqra8nbayUCgLOaa3jsvqHECK7k6QnRPZQA3Cz1RMbDuwGnkqm3VKgUWoX0qZeWAego1LqAOY5YTQwNIVTmgGLMhm3EDmK1NMTIhuwJadaWusrabSrBQzQWndzYN8/AIO11vvTbCxEDicjPSFczDaBRKeV8DCNtmAqk7s7qG8vzHo+SXgiV5CRnhBCiFxDRnpCCCFyDUl6Qgghcg1JekIIIXINSXpCCCFyDUl6Qgghcg1JekIIIXINSXpCCCFyDUl6Qgghco3/BzFrIkKw22apAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(7, 5)\n",
    "\n",
    "artists = brazil.plot_results(test_pois, results, ax=ax)\n",
    "invert_interval(test_pois, cls_obs, cls_exp)"
   ]
  }
 ],
 "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}