tutorials/3_posterior_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Posterior evaluation and analysis\n",
    "\n",
    "Once we have obtained a posterior estimate using SBI (see `2_inference`), it is important to evaluate the accuracy and reliability of this estimate before we can analyze it. \n",
    "The `sbi` package provides methods to do this: \n",
    "\n",
    "- we can sample parameter values from the posterior and compare the resulting simulated data with the real data to check the predicitive accuracy of the posterior\n",
    "- we can run simulation-based calibration (Talts et al. 2018) to check whether the posterior variances are well-calibrated (see main paper and Talts et al. 2018 for details)\n",
    "- we can analyze the posterior correlations to detect parameter interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sbi.analysis import pairplot\n",
    "from sbi.inference import SNPE\n",
    "from sbi.utils import posterior_nn\n",
    "from sbi.analysis import run_sbc, sbc_rank_plot\n",
    "from sbi.inference import prepare_for_sbi, simulate_for_sbi\n",
    "\n",
    "from consbi import DATA_PATH, RESULTS_PATH, BASE_PATH\n",
    "from consbi.simulators import RuleSimulator\n",
    "from consbi.simulators import default_rule\n",
    "\n",
    "plt.style.use(str(BASE_PATH.joinpath(\"figures/plotting_settings.mplstyle\")))\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load inference results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load inference results from previous 2_inference notebook.\n",
    "with open(RESULTS_PATH.joinpath(f\"tutorial_posterior_gaussian_05_N10000.p\"), \"rb\") as fh: \n",
    "    prior, posterior, theta, x = pickle.load(fh).values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior predictive checks\n",
    "\n",
    "To simulate data we need to set up the simulator first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the structural features of the model are saved to disk. specify the path:\n",
    "path_to_model = str(BASE_PATH.joinpath(\"data/structural_model\"))\n",
    "\n",
    "# set number of neuron pairs sampled from the connectome to mimick experimental settings, e.g., 50\n",
    "num_subsampling_pairs = 50\n",
    "\n",
    "# set the default rule as wiring rule (DSO)\n",
    "simulator = RuleSimulator(path_to_model, \n",
    "                          default_rule, \n",
    "                          verbose=True, \n",
    "                          num_subsampling_pairs=num_subsampling_pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "57aea53134634381ad632b9b12c14383",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Drawing 1000 posterior samples', max=1000.0, style=Progre…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "26686406514f4ed68d0563e17e6355f5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Running 1000 simulations in 100 batches.', style=Progress…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "num_simulations = 1000\n",
    "num_workers = 20  # number of CPUs\n",
    "xo = torch.tensor([[0.4300, 0.4300, 0.4200, 0.6400, 0.1700, 0.4400, 0.0900]])\n",
    "\n",
    "posterior.set_default_x(xo)\n",
    "\n",
    "# wrap the simulator to handle batches of parameters\n",
    "batch_simulator, prior = prepare_for_sbi(simulator, prior)\n",
    "\n",
    "# we simulate with parameters sampled from the posterio.\n",
    "theta_pred, x_pred = simulate_for_sbi(batch_simulator, posterior, num_simulations=num_simulations, \n",
    "                            num_workers=num_workers, \n",
    "                            simulation_batch_size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAJpCAYAAABM0i2yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAADOy0lEQVR4nOz9d3gk2X3fC39Opc4BOQ4waWdmZ2dzTtSSYg6iZcuiKMoWLSe99r02Ha7F59rPI/re97Upy9f3ubKvZNmWJQeakknJkklKJEWRy03cvLO7s7uzEzGDHBqdQ8Xz/lENDDBhZwZoAA3gfJ5nF4NCd9ep06eqvvWLQkqJQqFQKBQKxW5G2+oBKBQKhUKhUGw1ShApFAqFQqHY9ShBpFAoFAqFYtejBJFCoVAoFIpdjxJECoVCoVAodj1KECkUCoVCodj1GDf5epWjvxqxxvft2nnMVWz+848urNr29z50SM3jTfLc2QVeOLe4atsa53HXziHA7zx7nnzNXf5drcW1819+NMZCxQHUPK6VhYrNf1lxfVTn9LU5M1fmG69PX7H9kQNdPLi/a+WmG55DZSFSbCoNL9jqIewIbDWPLUHNY+touGou14tajzfOtdbbeuaw7QXRF77wBb7whS9s9TAULWIgHSUVvVnDpOJyDvYkt3oIO4KDvWoeW4Way/XTn46SjKjr442wrzuBoV1p/FnPOmx7QXT8+HGOHz++1cNQtAhNE3zyjkE64ibAVRe04vrs6YzzvkPdWEZ4CiuRuTYeu6WbA71JhFqG6+bRg90cVHO5LnRN8Kk7L10fFdcmETH42O0DywIyYmr82OEeBrOxNX+muMnWHZvum3ziiScAePLJJzd71zeC8pOvESklpbpHxNSImrqaxzXieAE1xyMdNdG0Nd2Kdv0cAlRtj0BKUlFTrcV1UrU9EhFDzeM6kFJSrLtk45Y6p69DEEhKDZdExMDUr2rjueE5VI+Vii1BCEFGPQWtG8vQsAxrq4ex7UkoN0XLUHO5foQQZOPqvL4RNK11c9X2K/e1iwVsz2fvF7+1vG0oG+PZL35gC0elUCgUCoViJ9H2gsj2fB7a38WTX/7E8raV4kihUCgUCoVivbS9ILoaQ9mYshgpFAqFQqFoGdtSEF0ufpTFSKFQKBQKxXpo+7R7hUKhUCgUio1GCSKFQqFQKBS7HiWIFAqFQqFQ7HqUIFIoFAqFQrHrUYJIoVAoFArFrmdbZpkpLtFwPOarDgOpKIaxWt/ank/D8QmbCy1VcxekowZCNRy6LkEgKdseCUsnkFB3PBCCuKVfq0Q8k4UaScsg06ycWrE9DE0QNfXNHHpbYXs+theQjoaVyR0voO76V12HDTd8LRKSUQPXD3D8K9+biW3vKudSSkoND5CYukbU0JgsNuiIWRga1zynFavPS+Ma52Gx7oIEXRckIwblhhvO82XnYRAETBYbdCUs4tbuuR1OF+poGkQMHYEgEdFpuAH5usNQJoqm7c51t3tWwA7kD1+b5A+PT1J3fBIRnZ+5f4QP39aPlJJnzizw5LtznJwuk6s6SCmxDI2hbIx7Rjv42LEB9nTGt/oQ2pYzcxV+cHKOiu0xW2rQcH0Wqw6aEBwZSPHjt/Zx/97O5defna/wa392mpliA10T3D6UYbQrzmzJRhOCQ31JPni075pCaicipeTJU/OcmCjiBZLuVITuhMXZ+QquL+mIm3z4tn4GszE8P+B778zy9OkFxnJVooZGKmqQjll0xC16UxE6V7y3M2Hxkdv66c9Et/owb5qLuRp/cmKaVy/mKTc84pbOfMXG0jQKdQc/gJ5UhGTE4HMPjfKBI71bPeS2YeV5GTV1HjnQxZ17sst/ny01+OM3p3lpbJFCzSUTM4iaOv3pKDHL4LbBNO8/3IumCZ47s8DvPDdGse5iGRofPdbP5x4c3bqD2wTGczX+7++d4q2pIsW6h6HD/u4EXhD2oIuaOp0Ji7/xvv3cPdKx1cPddNq+uWt05PawUvV7NHfd+8VvMbaikvUmsmUNDN+eLvJ/fONtgkBSsT0aro8mBJ+5fw8V2+PMXIX5sk2+5lKqO1QdHylDY1EqYnD3aAf/+OO3MtRxSRRJKXljosip2TKWoXHXniyjXYn1DvVGaKtGkBXb4zd/eJbxxRozxQa5qk3V9olZOg03nMeDvQkePtBNEMBCxeaF84toAhIRHRDMl21Gu+Ls704yXazj+AEP7e/i84/sveZTbQtoq0aQr48X+P7JueXfT8+WOTtfYW9Xgt50lM6ERdzS+dwDe/gX3z3F8+dy5KsOAI4foCHY153ggX2djOWqjOVq7O2K05uO0hG3SEUNfuHRfWhaS62dG7oWG67Pbz1znreniixUHAIpOT1bJmJo+IGkUPdAStIxk0zcRBeCz94/zMnZKjOlBkPZGJ+8Y4DHbulZ4zA3jZbPY8X2+O1nzuMFq1/y2QdG6M9ECQLJf3z2PCcmi0wXGwDNBxSIWQaDTfF8S1+KI/0p/uvzFyg1PGwvQBBe/44NZRjIxDg2lOah/V30pi8J7i26Pq7rnPYDyfHxPC+dzzNTbvDK2CJV26Nie5TqHhKJoQkCCRFDIx0z0TVB1ND5Sw+NYugaA9ko9+/t3M5WbtXcdafzozM5pIRCzaXqeEgZKvyvvjjOQCbKbKlBoe4SBBLb86k5ARCujLrj88zpBf790+f4/KP7KNc9YpbO+GKNl8YWWazaVB2fH51Z4LMPjXLncHZLj3WzOTNX4c2JInXXJ1exKTVcSnUPsyHw/AAJvDzm8NZkkbhl4AWSxapN1DJIRkx0DUr1sHt6vSmgAJ4+PU9vOspP3Dm4pce3WZyeqyz/ezJf49RchbrjMVe2ma/YdCUsYpbO3/ndHKdmyxTrDlUnCIV7833vzJTQhGSm5OD6PpqA8XydowNpBrMxppsioV3xA8n5hQp1J2Bvd5yZYgPHC1ishmIoV7Fx/ADXD6g74dqSwGLNpdxw0TWNX/v+WTRNIMSSqCzj+HLXWY7Oz1evEEMAr13MM9QRw3Z9yg2Xxaaodv2AcsNFAg03oFx3cfyAyUKdH7wzy4XFGpapoSGwPR/Hk8yXG2QTFi+OLTK2UOFnH9pLdzICwFOnF3j1Qv7SeBaqfOrOQQ70JDfl+NfCd9+a4anT87zZtNKOL9YIpMT1w/NMSrCRCMDzJXXHx2yK8//8/Hned6iPi4s1Xjy/yEdu62N/T3JHuxZ37pHtcGKWgR/IZTHkBQG+lHh+wJn5CjXbo+b4BHL1I5ds/q/U8Pjm65PMFBsc7E3heT6vjOcJfFiohjcfXdOYLdv8gw8f5thQZouOdPOZLzeouz4Qdp0OgvDi6niACC8iEknNDcjVPJYitOy6R7HuoTfv5jXbxRDQmYwghEAIwdm5MoWasys6WUeW418k08UGWnPuqrZLseFxbr5KxNAoNRz8QOL6l9bq0s+6G/DqxSIJS8PxJbYXkIgYvHwhzwOawGpjF2TV9vjay+Pkay4Auia4Yzg8j7xAMlNsUHc8bDfAv4pdxA3ADQIaXoAGaFqYBRMEkj96bXLXCaKIeeV3PZmvMVmoM5SNUXM8xhaqaM24NCnBDSSmFj7I1J1wLhcqNlKCL6HmBqyMsHQdSdVpMFtsUG44HOxN8vE7hmg4PsdXiKGlz3/lQr5tBVG54fLMmQVeOL+I7fpUGi51xye4ymsl4VwBOL6PAKYLNufmyti+xPFCcdmdjPDJOwbZ270pnoNNRwmibcqHb+vl66+MU26ElgjPlwRSkq9dejq6XAxdzkLV43tvz/L0qXkcX+IG4ZOCroVCwNAEfiD5zlszHOpLYe2SAM/uZISIoWF7AcmIwXy5QcAlMXk5l29aurn5Es4v1inUXTRNoydp8epFyT0jRX7scNu7PNbNXXuynJ2v4AehC8zxfAo1l9lSAylpurp0HE9eVRAsEQDlpoXTa3hIICMExZpDKtq+l7AXxxaXxRCE1qJ3pkt0xE1cLyBfc2hcQwxdTgAEzTtZrurw6sXFXSOsl9jfnSAdMynVwzm1PZ/pUoM7hrIAoeVCCCKGoO76CAGGEHhBeG2rO/5VLUyXbwnFQRjr9V+ev0jUMjg+XuSFczmycZP93cnla2Hd8TfwiNdHueFxfqGK7wdUbT98QL6B9y1ZKcu2x4tji3QmImTjJq4vcX3J90/O8Vce3bsjE3Pa92pyE1ze7HVp205u+KoJjcN9KS4u1qjaHgCGFi5kLwhwPHlDQQ2eBM+9dJqE7wdNhGqq5vi8djHPt96cxtAEezrj3D6UQW9t3EbbcHKmxOm5MsmIQcQImCnZobBcY6TNUlxI1NTI1wQIwTNn5/CCgJrjk4oa3L2ng0x8+2ZNVW2P4+MFclWHwUyUO4azWIbGns44n75riFcv5Hnx/CINN8Dx/WUBECzFzNwgglAY2K5PIhvjQG+S+bLdtskB04XGFdtmmrEtp+fK1JouspvFl1C1ff7n61P8pYdGd+SN6WoYusZfvG+Y58/mmCmFrsdb+9PL4sT2fKKGxkLFpuH61BwPx/dxvIBm4uJN4QWhi/LrL09wS1+KRMSgUHM5M1/h6EAagP097WspESJ8uC3UXYo3cZ4tEUgoNXw8v0G54S4Lz2Ldpeb4JCI7Qj6sYkcc0dWEz05v+Pr2VJGpYoM9HXHGchVsxycIQMoA+wbF0HsRyFAUaQIuLtZ4fTxPJmZxerbMKxfyfOBIL3s6YhsZILxpSCmZyNd54XyOsYUauibIxi1OTBZZrIYX1/XgNZ9OdSHoTUV4+lSOqXzoqgR4Z7rMzz4wsi1FUcP1+eqLFyk3wgvu2bkKZ+cr/PR9exAiDIoe6Yzz5LtzvDlZWLZyrAUJCAmWrtGVsLB0nWwbz1lnwmK2VKfc8HD90J391Ol5BFB11jERgC8lp2bKzJftVYG/O510NMxMhDCZ4b/86AIQlmM4MVlittTAC8J4oULNxVvnhXBlvNfBngSnZiuU6i6OF3BkIMWD+7rWe0gbQhBIKnWPs3MVSmsQQyupuwFeIDk+vsiD+7vpTka2c4D1e7IjBNFuZLJQp+H6BFKGMUQytEZcz012Mwghwswpy8AywgyrkzNlGq7PTLFObyrKn7t7iJ5UpEV73HyqtscfvDbJbLHOKxfyYYp8f4rupEUyYjBVrK/ZOrSSoJnhp2mCcsNt1qAJabg+r43neeLw9osJeWuqtCyGlpgqNLiQqy3HGbx0Pscbk8Vwfa5zf0vuDC8IuGc0SyravoLozuEMf3R8ksWqg+MFnJ2vIJDU1imwASoNj+++PYumCb70qaO7sm5MdzLCbYNp3poqMV+2qTselYZHxXYp2+ufYwnYvqTmhHGCMcvgzj0Z6m7Azz+6l/42FaLFussfvDLOU6cXmCs7674fhF4DyVTRZqpQ59N3De1YD4ESRNuMcsPlh6fmef5cjol8jVLdpeEFuDcSiHCTRAwNU9dw/YA3xvPMlR0cz0fXNfpKEeKWwQ9OzvHT9+9p+b43i+fO5lgo2yxWHaaLoWn41Yt5vCDA81snLiVQqHtUG2H8Q8zUGF+sMVe2kTLMBHz0YPe2q1NUbrhX3Z6rOJyaLfO9d2Z5eWyRQjUMnG4FUgYkLIPH2zz1/EyzxEAiYnBqpoyUkop9Y3Ec18OXMF+2+aPXJinVXd5/pJcnDvVuSyvjevjQ0T72dif4jSfPslh1wrgsrxUzHGJ7AWfmKsQsAyHA9yU96QjffH2Ku/ZkuW9FLbJ24cmTc3zzzWnenipdNWZqLWgCTE3Qk4rs6AQbJYi2EQ3X519//zSnZyuU6g4V26Paogvs1dAEaEJQsz3eqTh4fkDU1IkB08UGXckIFELL1HZ9YpjI13A8n7PzVfLVcE43QFsuc3ahyp6OGLYXUKzXl7cXai5PvjvPh472bdzON4DRrgSvXSys2iZEGCPz/Lkcr1zIU6y7eH4Yx9EKbE8ytlAjCIK2toyML9aRUuL6PguVBoFcv4VsJRIoNjzemizSnYzw0vlFfvzWPm4dSG/LgpU3wlJg+kyxQWfS4tb+FC+PLXIxV2Gu3MDxW381rLkBs6UGlq4hBHQkLMoNj6dPL2DoGof6kpyYLFGqu4x0xbmlN7mlcV3/8/Upzl2jRMFaCYLQTVtqrM/91u4oQbRNaLg+/+6ps/zw3fnlQoyVhrdhYgjA9UPzvh9IGs2fEklvKoHnS+ZLNkPZ2LYVQwCZmMlbk0VminUanr+hYgjCQM3pfJ1Sw6MrYZGKmfSlovSmIpycLvGBI73baj73dSe4ayTL6+MFpARDE9y/r5OvvTzOy80sK8fzcVu4UCUwUajxz//kJJ99YATL0OhNRdsuC1IT8MZkkXLDpVD38DbgZg0wW7Z59WIeIcDUNV6fKPCR2/q5tRn4u1OQUvKHr01ycbG2vO33Xhrn1HSZsVyVhhtm1m4EddslmoyQjpoU6y6uH2DqGq9eyPPy2OKy2/jNySK3DaaX45y2guliI7TGtnAuAsJaTm9OFPnvL4/z0/dtX6/Ae6EE0TbhrakiC2WbSjOjzPZuLF13PdTdACEgZoY1jzQRpl2WGi6lRljk7HMPjWzsIDaYB/Z18t23ZpYDXzeDRgCNqsti1WU4G2MgExYW3PQy8C3i/Yd7uWekg3zVoS8d5YXzC7w0lmex5tBwgg0R7b4v+YNXJ3lpLE9H3OT24QyfuH2Qw/2pDdjb2liKG6raHo7nt8xCdjlVx1uu85SJ1bg9lua5s7kdJ4jGF+urxNBS3aG5SgO/FYF+78FsxaXQ8OlNRehORpbjCifytSvi2N6aKnH/3k46EltTEqEjbjJbFNRbfEUxtLAP5h+/Oc2fu2uo7R5AWoESRG3M2EJ1ucrodLGB718KFLzJlitrpu4E+L6L2Vz8hibCdEsJtw6kmLpKavF2YrgjzoP7unju7MKGPV1eCwnMlhtMLtZIRQ1+7ND2sg5BWCX4RDNg+shAmguLVf6v755iqlDbMAEAYRLBYtWh7nqkIib5mkMQwEjnQWJWe2TA/Oh8DkHYQDMIgpa6MFYiAyjVXToTFqWaw/GmtS5qaDy4v3M5m3G7s1RjbYml7E9BWB5ko59nbC9grmwTt7RlMdCRsPCusuN8zdkyQfTE4R5mSg1K14jvWyumITB0jbrjU6g5OzK7cedJvB3CdLHOH742yfffmePiYo1CzWG6ZIMMzeKbZU6QhL5jKcOeN+mYSX86yn17OxnKxhnP11oWLLsVTBfrjOerm2YduhzHl0wV65i6tu0qD79yIc+T786zUHHI11y+8foUv/ZnZ3h3usQGeYdWIQkFe9l2mcjXOTFZ5EKuuvE7vkGKNTfsC2VqGxrrFBBaowAWaw6ThXA9zZVtvvnGdFvNyXoYvKxFSzJiNEsYbJz17XI8P1iunfXQ/i4ePXBl2r2hiS2N4bprT5ZMzGxJduxKijWXmWIYS9W5Q4P3lSBqU96cKFKou8stJOKWjq4JUlEDXRN46ynocpNoQmBoGtmEyUeO9nH3SMdyqn3c0tlmRo1lXjy/yO++OM6bk6UNi++4ERpuQL7qrHIHbAfemCgs/9tv9kkaWyjjrqEI3nrwfInrBVQdr+VPxethX3cCU9fw/IBGK4OoroIfSCq2S9XxsXSdfc2SB1LCGxPFDd33ZtGTivDgvk6W4pUtQ+fhA13LmZubgaELdCH49F2DPHygizuGs4ysKAyqCcH7DvVsab+vt6ZKxEyt5RbvhieZL9tETY1vnZjZNC/FZqJcZm2K60uCFStaCEF3MsJoV5yq7VGqu8xXnPf4hNaRjZtk4xZPHOomZq1+Mrhvb+e2rJRbd3xeOJcDIGbqLX+auhmWSuU/fXqeg73t2RfparhNETlVqDNRqDO5WKNQ8zZVDEnCc8XQAnJVh7PzVe5vkzX5k3cP8dLYIrmyu2HusiUCAAmZqMmjB7tW3ZDdLRT7reaRg93cOpBmuljnzFyFC7na5q43GVbE/m8vXOSe0Q4+dLSPv3DvMJOFOqW6y3BHbMtrY7m+ZDxf35DYvZilUW6EBR9X1hrbKShB1KYc7k/yznQJQxerfNSlusu7s+VNTX+sOT59aY3D/RnuHe3gxFSRQMLRgdS2jU/I15zlm1Q6ZoQ30C1SRRFDo1BzKNTCCrjbJVjxUF+KH56a5+x8WAZiptTY8DiOq2Fogphp0JeKMJmv8850maODWx9QLIFDvUkm8zU2w2ama4JHLhNDEH5P25kgkLx6Mc+ZuQoRU+OuPR2A4Ox8FS9YW/uTtbKUyXc+V+X8QpWpQp3PP7KXoWyMoctcelvFof4UuUrrYzt1ERaYnS3beIFkoWIrQaTYHA72pvixwz34MuCd6TJ+EBZfnC42mCvZ2JvlNAdsN2Cm2KA7abK3O7EjToLOhIWpC1xfUqx6GBotTQ2/UQwtdHd4gSQbNzk7X8HxAvb3JLb8SfN6PHqwmxfOLTJfsVmsOJsehyUASxfEIwapmMGtA2l0TYS9prZIEDnNQn6uH/D2ZJFSw92UGDtB2Nz0Q0f7GF+sM11sYBkad+/JclsbiMP18P2Tc7w5ecntdyFXIxUxAMnJ6fKmxKstIaWkZnscv5BnpDvOM2cW6EtH+dSdg5s3iOuQiRrUndavOb9pHbM9H5D0qaBqxWZy/95O/uGHj/Abn7uXv/zwXixdY6rZsuNmKTzzlTWPQyLxg4C3J8tr/ox2I2rqPH5LD0JA1XFp3EDTo/XM4dXQm4UvLUOjI26yULH59okZvn9yjt9+dozTs+0936au0Z+J0hm3EGHf2g1n5Xfw0L5Obh/OcqAnyZ7O+LKAjG6RhS1fdfid587znbfC7/DJU/O8MVG86tpq9Voym9mfyYjBzzwwwi/+2AH+5vv288jB7rZwH66Vhuvz9nRp1TYpYaJQJ19zeXem2BIL0Y1+H44PddenbHucma00i7pWKNbbJ3btyVNhv7zLacWaqzsBjhs2D2/Xpsrroa0E0aNf/j57v/itVf9FjPZIod0qdE0Qs3QmFusUamGQ9VouAMVnv7rmMQghqDgBX31lnP/03Bizpe2dar/Ens44iYge9iu7gdevZw4vx9QEyYhBdzLC+w/38MjBbowVmUh+IHny3flVcWTtyHI2zVIg1Aaz9B1EdXhgfxedCQtdE8t9pXRNcNee7MYP5Co8dzZHdUUPrcWqTdW5+sNLK9eSBvRnowx1xEjHwlTvmKXviMbLjh9cYWFz/YCFss0PT80zXbRbsp+b+T4CGVqKvEBSqLlhXFEL+tO1ij9+c/qq8UOtWHOWLuhMWMQtnfFtlgRyI7SVy2yyUGfsy59Yte2J5391i0bTXliGhi+DTXWVQWiKl1KGP5u+fMcL+IXH9m27mjmeH/D6RIHxxTqpqMGp2TJn5irUr3HT2igEkI6Z3DPawQdv7eV9h3r407dnr3hdxfYo2x6ZWPu6zm4bTHPXngynZssbWjX9cgY74gxlY/Slm0XyBKSiJveNdmxZfZSZVQ8KkrmyvSmCNmZpDGVj3DGcbZs4llaRjpr0pCLMly8Jn/MLVQo1h1LN2fSYtaVLXkD4oOp4knTMbKsG17om0IBWX9V0AamYSTZu4vqS2VJjx1mJ2koQKa7NSGeccsNDys2N/V3aVczSScdMcpWw39dEvsZo1/aKJfrmG9OcXwhrslRtj7enSzTczathsoSpCz58tI9//hfuWN7WlYiQuyxrMGbpJNqkyODVqNgeCxWbXMUlFTWpu615Wr8ehgh7qLVbU+HupEWp6ToZy9WotaDj+vXQBSQsg7hl7NjqwR+/fYBvvTHFQsUhkKFVZiJfp2xvfl8tjaXAaoEmBJ0Jk0/eMdBWbslU1AgH2uLrmqlrZKIGlqFj6mHW805DCaJtQK5i8ycnZpq9oti0G7ggvOBGTZ1M1FxVc2i7mePnSo1lMQSgaQI/kMyVGvibWNMJwnm93HDw4P5OLi7WLlXfFfDYwe62nec3J4p8/+QcE/ka08U6CUsnZgjqNxCLtV6GOmJhtfQ24+H9XUzk69Qcj9lSg3TMDFvsbOAJGzEEe7sTPLS/k5GunfW0vkRnwuIvPbyXfNVBCDg+XqBqe5v6YKiL0ErvBbJZfVyjPx3ln/7EbW0XXGxoGpYuqLfQOhk1BBFDAyEY6YyzvyfJ6A5cb+13VVFcwQ9PzVNuuPRnogRSMlVovGe4RuGZr1zVX3zhVz656vfMo58l+9jnrvoZS887S/uJR3Q0IehLR+lORRjcZt20L3+ajJk66ajBmC+vOpetmMOrYWhhXad0bPWp152M8JcfHuXkTAnbCzjYm6Q31Z5zXHd8nnx3jkBKHD/A0DUkYBgaeK2zilzrO3jql36cp4D/tzntv/zLv8yXvvSllu13rfSmo/z8I3t56fwiE4t1kFBtOMw++RUKG7CWQterxa0Dae4a6Vjv8NuejoSF4wV0xi00IZr+/Jv/nLWc25IwAcLUBUEgSUQM/sknbuXWwcwajmRj0YSg+OxXmf7Bf73ib2tZc5oIExgCKZatY7OlBi9fyHP/3s6Wjn2rUYJoGzCZr5ONm8yVbQYyMeZKNu57qP/sY5+7YpFf+JVPMvpL37zhfQoRngimrtERt5AS9nXHeeJwbzM7q31MxEs4XsCJqSJzJZueVIRjQ+nloPyhbAzL0JZbHEBYr2Ox5lCfCZtvrpzRVszh5YRp4mFm1sMHuq/4eyJicO9oe15g/EDy9lSJyUId2/Opuz6mrpGNmcwUQ4Fu6Rq6aF1PqWt9B5/5t8/x9z98mAf2td9cJSMGTxzu4Q+PTyKBvkwM74M/T6bFawkgYWm871APP3HXEPeMZNf1WdsFy9B4+EAXk4Ua9WmPmnPzdYhu9twWhNdCQxPhOrd0+jNRam0USL2S/kyEgx/9K0Qe+JlVc7PWNXegO04mEeFCrspgJkY2blFzfJ45vUDU0Ll9uP1E4VppT3u8YhXZhEVH3CITM/GCgGR0Y+NKNEAXAr3Zu+zHb+3jo8cG+MQdg3z02EBbuiv8QPL1Vyb44bvzvDNd4qlT83z9lYnllhxRU+fDR/uWYyyEgDuHs/zFe4cZ7owTMbSrpqq2CkH4ZJWKmjy4v2vb9S37xutTfO+dWd6ZLvH2VIkTk0U8PyAbt+hqNrg0dI3uVITEBrVziRjhhw51xHh7qn3bUZRtj66EhaVrpGMWg5koegvnI4wb0vn03cP88qdu497RjrZ8QNkoPnxbPx87NsCR/jQdCXNDWwdFdcjEDEw9bOiajZvs6YijaxrFWvuk2q/kn/7EMYazsZZczwwBuqYzka9TaYQxg9PF+vLf355u3/NwLbTfnU1xBbcPZvh/3pmj5njommAgE8PzqpSd1scmCMKgPMvQMHSNxw52LWc5iQ2VDOvj7HzlinIAcyWbM/MVjvSHhelu6Usx2pVYju/IxEzOzVf42suTWIZGICXuNVxo60EQPs1nYiYfPjbAX3t8f4v3sLFMFuqr4q+iZhhgP1e2GczGGO6Mk685+EHYbmaiUCcIgpbGEy0FUp9q2SduHOE5ZHLXSJaK7SEDyFWdlrTaWfISCRGu7xfPL/L+bSau10tnwuIXHtvHx+8YYLbQ4DeeOsPLY4sU615L4ystHeIRk864yVzFQdfCkAFNhDWfhjraM4YmETE41JfilYt57HWcgwLoz0ZIxwyODmZ4p1kPajJfpz8d3ZEiXAmibcDFfI0DPQnmK2EWT3cyQsTUeO1ia9W5BnQmLe4eyeJ6kpGuOB3xsK6JJgRHBtq3BcC1CqMVLnuKswxtVaroWK7K7cMZqrbLTMnGDySu799QocYbRQBeILFMva1jg67F1Z6ED/Qk6ExY1JyAYt0l0wwgzlcdNAHJqInh+TRsH7cFUyk0gdu823UnI9w21L5m+lTUZCgb5ZkzOfI1B0MLs5GKdRdnnf7EiC6wTI2eVIS66/HqxTwjXXEO9GyfHnitQIjwwXAgE+NzD+6l7oQlNSp1ryXp5hqXEh9cKcnEjHBNa4KRzgSD2WhbtIe5FjErjJHM19w1i8SIIUAKhBCkowbZmEmhHvbl8wKJqQtua8MYqvWgBNE24GKuSiJirHJVGZpGRBfYN3iBzTz62eu+JmZq3LUny7/+7D2UGi5PnVpgqlCnM2HxyIGutk6zvFb9leGO967LciFXI9OsCfSjszlqjh+WvRcBnidZqYtuZA6vhZQQt3R+oo1K/N8oQ9kYV7Z6C91/i9Vw/uKWzvmFKhKImgZdSZ1ExKBQc5jI19dV3kAjjN+oOj5HP/ELPLCvk/tG2zuIWBMCP7hUVLAjEeFC7lIhu7WupUzCpDMRIR01cH1JzfG4mKvtOkG0kkcPdvPOVInTcxUg7L14s21kLv8+IqZAF6F7XReCBw90kas66EJw90iWx27pJtmGoQNL3DaY4WuvTKzadjNrbumhJhk1SEXDXo8He5OM5Wo0PJ/ORPjgfKyNH0zWQvt+o4plkpFQ6S9RaXa7Nw0N1/dvqNzE9TIJoobgQF+ST905SNTUiZo6P3Xv8DpHvnkMZmPcO9rBqxfzSBm6FO7ak2X4OmbtVNSkUHPpSFj0Z2LMlcIAYTMIsLWAmn1pfm8mA2gJXYRPa8mIyf7uZNv3J7sambjJ47d088zpHEFTFR3pTzFVuBRLYOoah/pS9KQiZGMmQgi8IODlsTzlhsdwR5yLuSq1NTSMi5oaRwZSmLrO3X/7H5KKmm1trq/YHhOF+nLjYz+QvHxhkVTUxK2GxQTXspZMPazIHTXDy7bnBwQBJKO7+zKua4LDAynu3pMlV3WYLtSZLTduyl2Ufexzy+7IhKmzpyus+5aM6PSkIvSnowxkYhzpT/Gx2wc27FhaRUfcZCAdpVR30YQkuMk1l2mGFNw2mOFgT5KGF2aTHhvK8BN3De64AqBL7O4zaZvw4P4uvn1iZvn3iXyNgUyUquMTBHaz8Wvomnkvc7HGpQ4LglA0CAlRS6czYXHHYJaH9ndt6LFsJO871MMdwxnmyjY9yQgdCeu677l/bweT+TqGptGXjjBbahC3dIr1AMcL0DS4mTJFlh6KrJrto2mhUIhbBkMdMfZt46a49452cqgvxXSxQUfcoicV4d8/de6K12XjJp0Ji3zVxdA0DvQkEQL2dMZIGII3psrvmSF5OWEne41cxaEzEVal9ja5btTN4gdylTVNEv6ejYduM+nLZZEtAP0GaosZGmRjBh1xi4rtkas4GLrGmfkKp2bKHBvMEGvjIp4bjZRhjOCg7eEHElPXuJCrcj1NJIB0VF8+T01dYHtyOW4yYekMZmOIZsr9vW1umVxCInjoQBelhst0sUEQyOvOBYAODHXGeXBfBx1xi//PEwfIxi3GF+s4fsBoVxyzTWujtQIliLYBtw6kSVgGbzUza1w/wPMlVcfH88N2HlJKhICGG4R/b15gl4orIppuNlNDQ9DwfHqSEeIRnUzMoi8T5W8+caDtiozdLNm4RTZ+fSG0RFjxeJg3J4rc0pfk1oE0z5yeZ7ZkYzbr6whvtbvH0MIbdUBoWg7N84KOuMlwR5yBTBTH91moOFi6RmfSYjAT45GDV6babydSUXOVhetwf4pXLuRXveZIf5r3H+7l+HiBuXKDe0Y6OD6e562pEoZpMNoVZyxXvWI+UxEDSdg4E3lJtFuGhmUaeIGkanvELZ3Dfe0bywbh0/VgNspUIQzyN7SwdIUADCEYW6zhegGmoZGJ6tRciesH2E3rWUD48KJpoGsapha6Kz591xBDHVH+x6uTmLpGXzpKfybKXNnmqdPzfOS2/i075q3mUH+K4+MFEhGDkc44DdcnZhnU3asHWgtCq/hgR5yoeakifHcq0hT1EfZ0xMjEwjYVcUvnjuFsW7XoeC8O9Sd5+cIi+7oTVGyPqu1hIGhcQ3kLoCOmc+tghjv3dPDAvk7u3tNBJh6e7zu16OflKEG0TRjpii8vykDCqdkyh/tSmLrGTLGB4wf0pkK3z0S+znShTsPziVsGEUOj7gZhL5qoyVBHjM64RW86QsTQ6UpavP9w747rS3OjLAVnQti3rVBzOD1Xoe76CAQJS0cS3rQDKfElpKPhDTxhGWF2lZREzdDS9ksfO0JPKsL33pnj3HwFQxPcNpTZcUXMHjnQheMFvDNdIpBwsDfJB470EjV1Hj5wydI42hXj+HgBgJ50BFPXqLk+1YZHPGKQjhp85Fgfjh/w7OkFclUXz/PJxC0EgoodZld2Jy0ePtBFVxvHsi3x8dsH+NO3Z7m4WMPUNf7ifcO8Pl7ghfOL3BoxmMzX0bXQ6nDnnhQXcjUmCw10XZCJmuzpjJGJmpRsj95UhEcPdvPh2/qJmzpvT5evqNJ8br569YHsEoayMT50tI8fnc0x1BGj1HDDNh91l8WKjde0SuqaIGbqJKMGdw13MNodZ2yhRq4aJqzs607whQ8e2vZWkN5UlI8e60fXBPlaaCXSBDRcn6rj4/vhA7QQYdXtbDw8tz562wBPHO5p2wr5G40SRNuQRw50LdeC6ExYaBrcNpDmjckSmhDc0pvC9QMEkh873EvF9nh1xZO83uyq/hfvG6YnFVkuXqgIs1fuHe3glQv55fT7pWiVT945wGfuHyHwA4QQ1N2A333xYmjVIHyi//N371kWrj9x5yCuH6A1azrtNAxd44NH+3jicE8Yd3WNi+hoV4IH9nVRsV10IZgqNhhfrBHp0rhjMMNQR5yfum8YU9f4X98f4Po+3z4xw8XFOoEM3U0Sialp2yarJRU1+fP3DON4AboWfv+PHuymOzXOfNnm+MUCDc9nf3eC/kyMIwMZLEPw8w+Oohmh+2ZlEdGVPcqipn5FQ+KYuTtvYCs5NpThtsE0jh/w9lSJH5xcqqYueXuqSMP1OdibojNpcf9IluMTYRr5wd4k+4PQnf3wwa5tL4aWONKf5nBfis8+MMIfvDrJYrXBaxeL5KthCYG4paFpGulomFb/C4/upSPR/g8bG4kSRNuQjoTFzz+yl7PzFTxfsr8ngaVrTBQay2nmpq6RiIRm3obr8850edUFdjAbvW7A8W7lfYd6+YPXJlc1W41bOh+/fSAUj00BGbXgr71v/6rvIW6tPqV2ysX1vbje06RoZub86GwOCJ/mOxMWt/anuHukg9Gu+HKQtGWEBfAe2NfFRH6y+QHh/24dTLdlUdD34nIh83MPjjKWq7KnI8a5+SoR89LDyKMHe0iucPdeq1HrPSMdPHtmYdW2dq1wvtkIIYgYOkf607w0tkjV9olp4Zx5geTDR/s50JsgGTGYKtnMlULLkKYJIqa247KmhAibsP7Co3s5t1ClLx2uu9Oz5eUHvqGOOLcOpHe9GAIQ8uY65G1oO729X/wWY1/+xKptTzzxBABPPvnkuj9rA1jrY/+GzGPF9nhpbJG5UoOuRIT793Yu+4DnyzavXFikWHcZ7ohz72gHUbNtLENtNY8QNoP93ZfGuZCr0p+J8dP3DTPa1fZB0WuZx01rkXlissjJmTK6BrcPZZazsK7F+GKN18YLNByfA70J7trTsRmWtk1bi0tVv4WAo4Pp5QKiN8JbU0XemS6jidAycqj94qq2/Jwu1lxeGlskV7XpS0e5f2/nKkHdcH1eGltkqlAnG7e4b7SjHd2xLT2npZQcHy/w4tgi04U6vako9+/r5O492Z3sJrvhOVSCaH1s+Um/Q1Dz2BraWhBtE9RabA1qHluDOqfXzw3P4Y6VhAqFQqFQKBQ3ihJECoVCoVAodj1KECkUCoVCodj1KEGkUCgUCoVi16MEkUKhUCgUil2PEkQKhUKhUCh2PdurytlNMJSNsfeL37pi27Nf/MAWjUihUCgUCkW7smMF0dWEz+UCSaFQKBQKhQKUy0yhUCgUCoVi6yxEj375+0wW6qu2DWVjWzQahUKhUCgUu5ktE0SThfpmtNZQKBQKhUKhuC7KZaZQKBQKhWLXs2ODqq+GyjxTKBQKhUJxNXaVIFKZZwqFQqFQKK6GcpkpFAqFQqHY9Qgp5VaPQaFQKBQKhWJLURYihUKhUCgUux4liBQKhUKhUOx6lCBSKBQKhUKx61GCSKFQKBQKxa7nptLuhRAngMYGjWU7EpVSHtvqQSgUCoVCoVgfN1uHqCGlvG9DRrINEUK8vNVjUCgUCoVCsX6Uy0yhUCgUCsWuRwkihUKhUCgUu56bFUT/bkNGsX1R86FQKBQKxQ5AVapWKBQKhUKx61EuM4VCoVAoFLueGxJEQghdCPHPhRDTQoiKEOLrQoi+jR7cdkAI8ZtCiP+w1eNQKBQKhUKxdm407f5LwM8DfxnIAb8O/D7w2MYMq/0RQgjgnwJ/A/itm3y78lOuRqzxfbt2Hi/mavz+qxOrtv29Dx1ayzzu2jkE+N7bs7w5WVz+fY1zCLt8HgG+f3KW18fDuVTzuDb8QPJbz5yjavuAOqffi4rt8R+fOY8frD7cn75/D0PZ2MpNNzyH1xVEQggL+LvA35FS/mlz288A54UQj0gpn7vRne0UhBD7CUXQMeDiFg9HsQs5O1/Z6iHsCNQ8to6zc9WtHsK2Z6bUWBZDivdmbKF6hRgCODtXuVwQ3TA34jK7C0gBTy5tkFKOAWPA42va6/bnYeAccDtwfovHotiFRE19q4ewI1Dz2DqipgpJXS9RQ83hjXKtc3c95/SNzP5w8+fkZdungD1r3vM2Rkr5FSnlX5VSzmz1WG6EL3zhC3zhC1/Y6mEoWshtQ2ki6ga0bu4Z6djqIewY7lZzuW66khH2dse3ehjbgn3dCToT1qptUVPn6GB6zZ95IzFEcSCQUrqXbbeB6Jr3rNg0jh8/vtVDULSYdNTkp+/bw0vnF1msOQyu0US827l9OINpCN6cKKIqkKyPY0MZLEPjjYni9V+suCafuH2Ql8cWGcvVtnoobY2uCX7q3mFePL/IVLFOV8Li/r2dJCM325HsEjfyzjqgCSEMKaW3YnsEUE5jhWKL6E5G+NjtA1s9jG3Pkf40R/rX/lSpuMShvhSH+lJbPYxtjWVoPHKwm0cObvVI2p9ExOD9R3pb9nk3YnMfb/68/Mo7yJVuNIVCoVAoFIptx41YiF4HysCPAf8VQAixF9gLPLVRA1OsjUe//H0mC/VV22bO5YgYKnhUoVAoFIprcV1BJKW0hRC/DvxLIcQCMEdYh+iHUsrnN3qAiptjslBn7MufWLXtied/lefP5bZoRAqFQqFQtD83Gn30TwCT0EJkAt8G/vZGDUqhUCgUCoViM7khQdQMpv4Hzf8UK5BSPrHVY1AoFAqFQrE+VCEThUKhUCgUux4liBQKhUKhUOx6lCBSKBQKhUKx61GCSKFQKBQKxa5HCSKFQqFQKBS7nrU3/VAoFKtYrDqcnCkxX7bpS0c43Jem47Lmg4q1UWq4nJ6tYOqCQ30poqaOH0jOzVco1F2GsrFt38+t7vicnivjBZKDvUkihsbp2Qq257O/O6nW0k3QcH1OzZZx/XAuMzHzPV/vB5Kz8xWKdZfhjhgDme29lq6GlJKLizXOzVWpui4jHQkO9afW1R1+p7FlgkgI8QTwA+BjUspvCyE+D/w2cKuU8uQG7/tXgcellA9t5H4UOwspJRXbI2bqGPpq4+qJySLfemOaE1NFHC/AMjSODaX55B2D3DaY2aIRbw88P6Du+iQjBkKIK/5+dr7Ct96Yxg/C7qvPnsnx6bsG+eGpeWaKjeXX3TvawfsO9WzauFvJXLnB778yScP1CaTkj9+YxtAEkebN6pnTOT58Wx+3Dqiea9djoWLz9VcmqDs+AM+eWeATdwxwoCd51dc7XsDXX5lgtnRpLd05nOXhA13ErJ0hFqSUfPONaV44l+PMXAUJdMQt7h7J8NP3jVwhthuuj5TsmOO/UdrJQvQt4GFgbCN3IoT4R8A/BF7YyP0odhYT+Rrfe3uWfM0lYmo8uK+Te0c7AXD9gKdOzzOer+F4ARBeZCcW6zx9eoHDfakrBJQi5LWLeZ4/t0jD9cnETD54ax8jXfHlv0sp+eG788tiCMKL9ddeGScIVn/Wqxfz3D6U2ZaWlOfO5Gi4PnOlBhfzNeZKNoGUPHqwi454hEBKnjo1z6G+FLp2pWhUXOLZMwvLYghC68+T786zvztxVcH95mRxWQz5geTcQoUXzud44XyOw/0pPnqsn7jVTrfKm+dCrsbp2TJjuSpLZ1K+5jBZaPD8udxyk2jb8/nTt2c5M1cBYF93gg8f7d81wqhtvmUp5Twwv1GfL4Q4BPxfwMeA4kbtR7HzcLyA//n6FLYb3oGLNZfffmaMP3h1klTUYLQrzmLVoWp7q95XdXzqjk+p4dG5DW/SG81Evsa3T8wwma8zW2oQSMmrF/L80scOs6czAYDtBRTr7hXvHc/VGepY7daQEubK9rYSRItVh1cu5PnTt2eQwEyxQcX2WKw6ALwzVea+vSaWoVFzfMoNl2x8+xzfVjBftpf/LaVktmzzznSJzrjFA/s7GbrMtTpfblB3faYKdSbzNRpuQDJqUHM8LuRq/Onbs3z6rqHNPoyWEASS4xMFvtcUORXbw/UldcdHEzBXslZZxp4+tcDp2cry7+fmq/zg3Tk+fvvlvd13JpsiiEQoy/834BeBAeA54D9d9prPs8JlJoT4EvB54G8B/ww4ApwnrJZ9Bvg3wGOEIupXpZT/5jrD+HdABngC+P8C0XUfmGJXcHGxtiyGbM/nxGSRM81Yj75UhP5sDMeXxIzVVqCEpRM1NRYrNlOFOqNdcVLR945l2Mk0XJ9z8xXmyzaZuMmZuQonpopMLtapOh66JlisOvz7p8/zv7z/IN3JCJP5GjXHw9Q1zBVWtuGOGPIq++hJRTbvgNZJse7yuy9dxHYDPD+0TMyWbKKmRiAlfiCZLNTpmK9wsDdJKmqwULGZyF+5llw/4PxClUBK9nUndnUz555UhHIjfDh57WKeqUKDTNzk/EKFd2ZKPLSvk6ODGfoz4S0gZhq8OVGg5vjMlhoIQiG+ZBU6v1DF9YNV62+78N23Z3lnukTF9ijUHCbzdTQBgQzF0rszJQYy4XnmB/DMmQWiprZq/ZyerSCPyata13Yam2Uh+mfAPwJ+BXgS+Ajw72/gfb3AvwX+MaHw+ZfAV5v//k3gXwF/F/jXQogXpJQvvcdn/R0p5RvArvhiFa3D0jWklAQSphbrvN68eGoCyrZPrupw/75OPF9i6gLXl1i6Rn86Srnh8T9fn1pecx862sexod0XUzRVqPP7r07w+niBcsPD0jVcP+BCrkrZ9hCAJLyxl+ouz51doNzwmC/b1Byf8wslDvQk6ExEiJgan75rkCffnWeudMkacNdIdltZ4k5MFqk3rYo9qQgvjS1Sdzz8QEMIgZRhMPn5hTDY90BPgm+8Pg2AJgQfOtrH0cE0uYrNH7w6SaX5WRFT4yfvHtqRgcE3wsMHupjM1/jO27Oca8bLlG2P7749QyZmMTZf4chAmlsHUnz02ACO71OouRTrLnXHJ5Bg6hpeEAA6pq6hbcN7RrnhcnKmBEA2ZpKOGviBpOz64fUskJRtj++9Pcvb02U6ExblhovrS/Z2xelvrh/TELvmnrnhgkgIkQb+PvDrUsr/vbn5u0KIFPDXr/P2GPB3pZT/o/lZfcDvAL8mpfzV5ra3gIvAI8A1BdGSGFIobpYfnVvgW29Ok6uEN2fXlyxdH4IgYKHq8u50hXtHO/jMfXtYrDl0xCP86FyOVy/k8QKJF0gsXfDGRIG/96FD3DGc3dJj2my+f3KOi7na8pO77flM5Gvkqg5SSiThU2vV9vmTN6d49myOvlSEiKFj6IJExMB2A953qIdbB1LELYPP3LeHM/MVCrUwM2i4I/7eg2gjao7HV56/wJuTRfxA4ng+VdvDC8B3A0xNoGkCXQhcL8D2fMZyNbqTEYQQBFLyg3fnONib5OnTC8tiCMB2A35wcp6ffXBkC49w6+hNRTnYm+JrL08QSBBAw/UYW3Ax9Tp+IHn5Qp6+TJR83eXl84ssVGzqbiiGdC0U5lXbIxkxOTaU2ZZxW3UnDIz2fZ8nTy0wtlCh5gZoAnQh8AIQInQp1lyfuGVg6gIJeEFAVzKCqWvctYuuVZthIXoYsIA/umz7V7m+IAJ4dsW/Z5s/n1+xbaH5s2NNo1Mo3oPfe/Ei//7pcxSaMR1+EN68pQQkTbeNZK4SBidezNf4Sw+OMFOqcyFXxZdyOT4mHTMRQvBHx6foTUWXTfY7HdcPmC/bzJYazJYalBsenh/g+hJDE9heaH1bcoH5EvJVm3zVJmYa7OkMn1SFgNGu+LIrw9A1jvS3d9ZVrmLzxkSRuusTMUKrmKFpPHd2gfF8DT+QVG2Piu1dmgMJji8RvqQnZRExdaYKdXRN49aBFMlI6CpzvICFis1UsX7FfmdLDTw/2HXB/EEgeXu6xB++Nkmp4aKLcE4dT9LwAmwvFDxeIHFyNf77i+PYnk/d9cNzu3lel22PxarLIwe6efxg91Yf1k0zXazz5kSRiXyN1y8WmCjUWTrD/AA8JILwIcTxJfmaiyYEmtDxJUwXw2vW+w51c8/I7rm1boYg6mr+vDxgevoG31+6yrbqin9fLZRAcRkRQ2fvF7+1attQNsazX/zAFo2oPXH9gIl8naipcX6+ym8/e55c1cbxQquQLsIbNsDKJCfb8ZjzAhYqDf5FoU7MDGNeclWHmuOja1C1BZmoScTQeHOyQM1Jkoqa2yruZS3kKjan50q8O1MiV3EQzRgGAU0ryKU5BfBWTGzN9jgzX2FfZxzH85lYrNGVsLaFCX+21OBrL4/j+pLpYp0LuRqD2SgjnQm+89YshpA4foDj+6uOX674mSs7SKBqewRScvxinntHO6k6HuWGx5nZCh0xkxnXXrXvTMzcNWKo7vhMFetkYyavXMjz1lSJXMXG8YJLIruZkhhaP8J/Ob7P+bkiaBqOF+D5l87puKVhaHBypoyhz3DPSMe2eYA5O1/hG69PISXkKk6YWSbDY1t5s1z570CGAf4VWyMd0TEMHccL2NuV4PxClWTEoDe9PY5/PWyGIFoSQv3A6yu2bz/ZvY25eyTLk1/+xKptlwuk3c5EvsY3Xp+m4fpM5mu8cjHPfNnG9oLltO/gGvLb9sGQAQGC8XwNGVy6uErA0MJAzUN9SYp1l++cmGGo6eI50JvkE7cPbEuz/HuxlJ33rTemeOH8YmjCh+UrsQZ4/ns/z/hA4EnOzFcZ6ojx7bdmeHOyyE/eM0wy0jZJslfl5bE8ri8JZBgcDWEWWanuMleqhy4L3vuJzgPmyqF1UhPw9nSZ8Xwd2wstTSeny/RnIhzqT5NoWs6EgMdu2R2X1xOTRX5wcg4vkNiuz3i+zoGeBIMdMc7MV3E8H0PXmkkRV8500YHVjzYhdSdYFlcvX1jk3ZkSwx1xPn3XEJbR3kLzhXOLSCk5PRsGkXs3aDLwJdTdgIYbkIlJCnWXX/n2yeVYtNGuOJ+6c3BbBpffKJtxRXkOqAGfAb6zYvunN2HfCsV1eeXCIq9dLPCjczkSlk5PMsqFxRqLVWc5mPJaQmgJCbjB0r+uxAug0vA5NVshV3G4b2/n8t/OzlU4MVnkzj3ZlhxPu/DqxTwnJgu8O1tefkJfyZVbro4kvFhfXKzz2tgikYM9PHN6gY8e62/peFtNqeEyW2owka9xMVcjaupoAk5MNvCC8LhuxrwtmzesuhsKJF2EsUjFhosXwOcf2UvM0jncl9oVT/NV2+P7J+eWH1bspgsxHTWYL9v4QdCMzfJveK0tIQmFaCoaWkoius7LF/LMlmz++vv2tXUWX7nhkq+5TBXr1Gz/+m+4DAkUGx4vnJ1nsCNBXzqKJgQXcjWOjxe4f8W1a6ex4YJISlkVQvwy8C+EECXCAoyPA39zo/etUFyP4+MFnjq1QMP1qTQ8Kg2PqXyd2bJN3fERAgwtzBxbLwGhGyV+lSJnFxdrO04QXVysMV1ogIQgaI1v+4ULBTqTEVLR9rYOAUgk5xdC776uCfI1h1Ldwbn5e1Tz81bjS2i4AabmMZGvMVWo8/OP7N0W7sRWMFWoryrYmYgY6JrgxFSRuVID1w9Ygx5YRhI+xPiBpOH5aELw/LkFKrbL+w71cNtAhky8/cpo7OmM8850CSnD5IW1ICXkqi4BdZ58d45DvaHIvpirKUG0XqSU/1IIUSbMNvtFQtfZXwd+dzP2r1BcizcnClRsD8fzcbyAuutxMVdbzgwLVvp4WkAQQL7qMlGokY6YJCIGqahJ+jq9lrYDnhfw9kyJjoRFKmKQq9iYRnhz9q9nYrsJnjw1x96eBBdyVeKW0bYxWBqCVNSg3PCIGBozRXfNYuhaSKDhBURNnXzNZSJfZ0/n9sm2Ww9L50wgJTXHwzJ0DvQk+faJaQo1F7cFay4gdOsKQotf3fGx3QBNCF67WOCn7h2mLx2l3HCp2B49yciWx27d0pvkpfMmL1fsdZ13DU+Sr9pUGi7FmkM2btGVMLiQq5KJmTuyQOimPWZJKX+TsHbQSn5vxd9/hzClfun3LwFfuuwzvk3odl+5rXH5thsYyxM38/qdylA2dtU4ot0SbD2Rr/HU6XnOzlWXL3aNZurtRuEHklLD4fmzi01LkWCkM8bnHxnduJ1uAq9cyPPrPzhDxfaYKzewvYDeVITFikOuat+0y+K9qLuS333xIj86k2Nvd4IH9nXxE3cNtl1MkaYJbhvMUG64vD5eQBfXixhaG4GUHGz26ZK7KMWkLx0lbuk8dWoeLwhFSyyiU6g5Tfd1a3ADyVypQb7uEjN1pksNpgsNRrri/OjsAqmoyYnJEoGUxCydDx/tY/81+qZtBt98YxrbC8LegOv8rIYnsT2fquMzVWxwbqHKN9+Y4fBAmkcOdPHx2wd2VExRe11BFJvKtUTPdg22Xqw6/OhsjtlSg950hIf3d9GVvLr1wA8kf/zmNHXHp+Z42K6PFwTLN5SNuq8EQMOVCOETNTVSEQPHD3jq9EJbtwcoN1yeO5tjqlCnI27x4P7O5WBLzwv4jSdDMVS1PWZLDVxPkq86CCmxbzSq8yYo1T3O56rkaw4N1ycV1fnUne01f7f2p5nM14lZYTPgVorClfQkIwgBqahxRTuTnUzD9bE9nz2dcUp1F0MTvHIxvypLsRUEEuYqDhqSmgRTDzg+nkcIsF1/uQEvhBlvf3Jihr/++P4tDb6uOj7I4IYDqt+LpTIjmgDXl4zn60RMnUzUoCcZ4ZFtWJbgWuwcaafY1TRcn6+9PM6p2TLFusvp2QpfW9Hx+nLmyg2qto/nBwSBXK5Bomvhib+RSEIzfBCA4wcYmuD18QJ1x+f4eIEXzy+Sq9jX/ZzNwg8kX39lgrenShRqLucXqvz+KxMUas1+WzMlZkth9lSuYtNwArxAUncDahsghiCMn6k0wtTz4+MFvvLCRU5MFpeb67YDtw9nuHM4w/n5KotVh40I7dEI6zHNV2xuHUgR7CIT0USz3UTMDIt3ThRqzJYaG/IwI6XEDSSOFzRr9wS8M11ioeJc8VrHC5azCrcKXYPxfGvHsCSKHM/H9QMKdZdzC9Xrv3EboQSRYkdwerZC7TLxU3d8Ts2Wr/r6uGUwXQyDpyWXRFA8YqBtQvq760sqtke54WK7YZ+k//yjMX5wco5nzyzwX56/wNtTVyvBtflcyFUp1FY3WHV9yVtTJTw/4EdncyxWXUoNl0LdvaLeyUbhS5ZbLpycKfPvnjrHf3vhAjXHu/6bN4HxXI2vvHCRhYqNIKDmtF6sBYRFLDUhePF8nq++eJGG2+JApTYlbhmM5aq8djHPy2OLvDFRbLl1aAnbC9BoVrH2wgrQMUsnYl79FpqIbF0WWsX2eGeqhOO29iwMZNgyRgiBJgSmrrWdm3q9KEGk2BE4/tVvAo5/9StkxNCoOz6pqImpC3RNQ9cEQm68hWgJLwhwmhaU/nR0laCTEp4+Pd/SYOS1cq0MO8cPODlTpur49C+neW/ueJe+XUvXmC01mCo0OH6xsKljuBa/+/I4Fduj7vrMFDfO4tfwAmTTMpSrOLw+XtiwfbUTMVOnanssVh0WKs6GiSG4VNHZD0A2g+Uf2t/FSGfsCtfYaFec3tTWlT24sFBloWLTalksCS3xCUsnEdHpTlrcO7qzqlgrQaTYERzsSV3RgFETYjnY9HIWqw57uxMc6U+xvydFbzqCoUHF8a4Z81J45istHXMgIRnV+cUnDlzVKlVz/FU9qraK0a74VeMhbulNMt907T1yoJMj/Sl0sbGXlMu/A0HYlynTzDiqOt7ymLaCc/MV/vC1Sb728jgnJotUbJeLizXKl+V/t3ItBXIpPTxUBFdz4+xETs6UMLSwGryzzrIYN/J9LLnKHD8IK0AvVLhzTwefuX8Ptw2m2dMZ59GD3XzqzsF1jWW9FOpumCF7+fYWrDlTD8Xgp+8a4qfvH9lxGY1KECl2BJm4ycdv71+uT5OMGHz0WD8d1+h+3pmwiBgaA5kYh/sT5Koutiff8ymz+OxXWzZeSXgjNzWNbMyi4yr1TBIRvS1M0lFT5yfuHCTbHGPM0nn/kV6GO+J0xi0WKjYLVQc/kPgbHMOy8jsQhC6MpVpRENai2aqn89OzZf7o+BTnF6pM5OsUqjbji3WqtneF3ayVa0kQuklOzoTu4URE542JAmfmKm1hYdwI5ss2z5xeWA6qXy838n2s7GFYbri8MVmi0vDoTkb48G39/NS9wzywr3PLs64GMtGrVoBvxZqLWQZCCJ441MtQducF8G/91VahaBG39KU40JMMOzeb+nvGAkVNnUcPdvPku/O8OVGkWHNoQe3Fm0KIMBvkD16boDcVRUq5XMNEE4IfO9TbNu089nTG+fwje6k6PjFTR9cECxWb584uMFe2KdZsTkyWcK/homw1grAdSiZmYhkaVcfncH+K4Y4Yd49kN2UMl/PKhfyq35Mxk4brs9FTslSAb6pQZygb45UL+WVraXfS4qfu3UPsKsVAtzOvjxcwdI3ZYn3Tz1tTF2iaIKJrfOftGT54tG9zB3Ad9nTGr5lMsh4EYXB53fU3JEGgHVCCSLGj0DRxw1aVu0c6GOqI8a03pjdfDBE+cUZNjYYboGuCbNzi/r2d2F7AwZ5k21XBFWL13D5zeoGGG2BoMFVsUHeDTY0g6ktHuW9vB14Qxn394hMHONiT3LLCeJcH9VcbHqmoQaXhYW/gAouYgs6EhSYEgZSrXMcLFYdXLuR3XG+zuuvzo7MLjG9yNlfYHFYS1TUMQ1BtA5f25cSssEVMq09GSRhPOJiJXTOYfLujBJFiV3N+vrpcTXkzCS8uYcXboWzo4inUXA73p7bc5H6jTBcbzJcbvDERZsMJsXmFAQ1NMNIVZ09nAk0IDvWlONKf3pydX4P9PQleWxHQrWmCqKnjB2DX3Wu/cZ0Yuk4qauIFwVWLik4XtzYFvJXYns/r40WefHeOE1OlLSlEKSVEDUHCMjg2lNn8AVwHU9fWXZDxWkgkx4bSV8Rr7hSUIFLsak7NlslEzSu6jhee+cpVfe4XfuWTq37PPPpZso99bk37DoKwx5XVbBSZiZnLsTDbga6kxTOn51maOU1cvwnuzXC97+AM8FXgp/7aF/gb/+ZXW7fjNfLwgS6KdZdz82FtluGOGMWaw7nv/A4zT/7XK17fqrWUiuiYuuCW3vRycPlKuq9RnHS7EQSS339lktlSg7NzFQxtbQJ8ved2ICWarnPbYIa//NDemx/ABjNbrDP31FcobMD1KwhgtmTztVcm+NkHRtrGpd8qlCBS7Go0TRCP6GED1xV38+xjn7viQnHhVz7J6C99s2X7jpoC1w+YKzfoT8d47JbubdWY85EDXfzxG1PLv+uawGuhInqv7+Dh/Z38xs/dG76uTXoqRQydT981RKnh4vmSb74+xdtTJYY/8HNEHvyZVa9t1VqydMFANsbRgTQPH+gmZun86Gxu+e+JiM49OyQ1eixXZbbUAEDXNFJRE8cLaNxk8c+1nttaM3g/GTH44JFefuaBPSTbsMnwRKFO9/s+R6bF1y9dQCpmkKs6zBUbnJmrcLg/td7hthXt920qFJvIHUNZfvDOLJYuWtIM8kbRCa0pnh8WZfzcQyPb7kl+uCPOP/rIEX7p99/A9gLywmEz6hCZWhgU3y5C6HLS0dBKcz5XJRu3wo7j1Y2JNTF1wcVcjaFMGEyejVsMd8Q4vxA2vj06kN4xAdUV26Nie5ydr3AxVyVfdxGbGLUmCL9bIULL8n9+bozHD/XwxOHetsgGXcLzJZZOS/u5AUQMQVfCwm82vq7YG+cG3iq2R7CCQrFB3D6c4XB/Gm0TzwQByGa8janry0UityNHhzJ86dPHwgDwTbRubYdkck0IBKHo3QjipoYmBA3X5/WJIv/jtUmCQDLcEefxW3q4d7Rjx4ghCIuXnpwp8e5MmYbnY2hseAbfSpaC1pesoMmIyenZCt98feo679xcDvUl8YLWnouWLuhORsnELOKWjmVojHQmWrqPdkAJIsWuZqpQJ19zNrVei2z+TwhBbzpC3DJ4q03adKyFVETnjuEMB3oSxDYj+0SwYUGjreT2wTRuEFCobkyhRF+GPfgCCWXbpVBzt7yH1kayUHHoT0eXM7t0TWxatpMuwp5xjh/Qk4zQnYwsu8umi4226j34C4/tb+m8hOUtjOWq/rf0pXjslm56UtvLon0jtI+dT6HYAl44nyNfc5FXhFVfSebRz7ZsvwlLY09XAlPXmhf47WDzWE3V9pqFCCu8MpYnV3E2/DAyj36WmGm0fTDndLHO+VyVuuNfdUpasZbCiuoSXYRPtotVe0uyrjYLiaQ/HSUVNcLmyFK2rG/be30fyYhObyrKgZ44k4U6ccskV3Vwp0oc6gvLPLTTtCcjBrf0pXjlQmHV9rWuOT8IXdSH+pMMZmP89H17GO3aedYhUBYixS6nUAs7tMsbsBCtNZvschKWoCcdJdWMNbG9gKMD7Ze+ez2ePj3PbKlBzNSp2GGxtnQ8fIrcKLrf9zlipoYXyLZ1M0op+dYb01zI1TA0gbiKP7ZVawlCT6Vl6kzk63Qn2zOuaj24fsBTp+Z55vQCp+YqpGMGhh42F21V+ve1vg9dwP7uBB871o8XwP6e1LIYLzVcJgp1+tLRtor/e/r0AnDlzX09a64/HSVm6hiaxjNnFtYxuvZGCSLFruJirsYzpxd4c6KI4wUMZmP4gdzQwnkr0QXs7U7Sn4khCBtUfvquQUa6tl9PoLFcDQiDXU1dICU4XkBmgzJvBGFpgo5EhLip80fHJzdkP+ulUHN55UKeC7kadde/ZnPcVmHqGh1xk1t6U8w0s7B2Et99a5ZXLuSpOT4jHXH6UjGycYuIqdGXiRLZQAEes3T2dSeJWzo9qQgHepLsbfb2EyKc+0/dObBh+18LF3JVOuJmy5pUL7XHWUoWmCvZ1Jz2K0jZCpTLTLFr+MG7c6s6ob96Mc8D+zqYLNY3zeQ93BFjpDPBLb1JdE3j/n0d3NK7PVNXkxGDajPrp9RwQ0uFrhE1NXJVt6VzKoC4qbOnM04mZjHSGWe62GCu3NjSzuJXI5CS8XyNfNWhZnvIDVpdpiboTJp0JiL0paPErPbofddK6o7P6bny8u/JqMGde7J0pyxcT1Ksu/z+KxNMb0ALD0NAJmpyS1+Sn3lghP/4zHkA+jMx+jNhH689nfFlS2+7kIwa2G6AZWh4LUg1swydka4kEVNv/q5hbZPisTfLzjp7FIprUKg5vD5eWLVtserwJydm6U1EyJVtvGBjQ2AiOvz3X3yEQErmyzZvTZV44dwiL5xbZCAT5RN3DLTdxfW9eGBfJ7/z7BgNNyBhGVQdj86EScX20QQtu0EJIBXRiFo6e7sSjHTGl+s1Xa2J5VYjYTlIX7IxAeC6gIipYek6CctgKBtjX3eC3nR7icP14gXBVeOixhZqy+KvM2lRqDvUbP+KDu9rRRDWKNvTGeOd6TL/6bkxposN8jWHW3pTzfYYgvv3tl+Np/v3dvInb0yTiBjUXWfN17SkpWMaGt1Ja1Uj13tGOrasPc5GowSRYleQqzpXvbBOF+ocGUxzPlej4fp4QWuf5wVhQTdNCEa7EvQ1b1hvT5WYKV5yb0wXG3zvnVl+8u7hFu59YznUl+KukQyFmkNH3EQIgZQSN5DUHY+FamvqlERNEJqGqev0pqPLYigdM+lvQwGgCcFAOoouBH4gKTW8GwjZv8HPBmKmxkhXnAf3d6JpGoPNGkR3DmdbsIf2IhU1GcrGrsieS0UNZLOOV8zUiZo6Nbt1MWURQ7CvK0EqahFIiRCCPZ1xYpZOseFyz2iWO4azDLZhx/dDfSk+/+g+/n9//A6Fmrum+moRHfZ0xbmlN0VX0uJwXwpPSg71Jbe8Rc5GogSRYlfQm4os1xFZyWhXHNeX3DOS5eULeTzHRxD6zb11PG5mojqOH6ZFCxG6Nv7m+/Yv//3sfOWK91zI1fD8YNs8fc2XbQxNoytp0RG3loVKw/UZW6hSang467DgWHoosEAjHTW5bTDDuzMlbhvMsKczzoeO9qG1YbZZOmpwx54sL40tMlu2iRgani/XVfhTA5JRnYFMjCP9ae4ayfJXHt3XukG3MR+7vZ/vvjXLxcUapi44NpQhYmj86GyOd2bKBFKSjhpUGu5NV62+HFODnlSU4Y4Ye7sTSKB3RcD0UvD0Iwe7l2Nq2pEPHu2jYnv8qz99l4uLN1eKIWlp7O1Okoqa9KUj7O1O8PE72itOaqNQgkixK0hFTR7a38lzK9oa9GeifPDWXv7wtSluG8qQq9iM5+uYhsDSdTQkvZkoh/tSLFZdjo8XqDreciq1IcJqsEuXYE2ErgxTgw8f7WekK0Gh5iCEoC8d5QO39i3vO2rqVC97orWM1mXNbDQvnl/k2TMLBFIyV7KZLNQ5OpBG1zTuGMrw4vkciYgBtrdcKyes5SIQUuL4EBDOnU7YFy0A6o5HJmoQtUzipkau6tCbjvLwgS4iho6UkocOdPHIgfbt3m7oGh87NsDJ6TK9KQtTE0gkVdun3PCW3WiXf9Mr9XdEFwxmo6SjYYr3/p4Etw9l0LRwjcR3UMHF65GKmvyFe4dpuD6GJjB0jYbr88rFPFXbQxNhP8ChjhjzZRvHD0sRhBWbw2y0UsMjkFzhUktGdIJAhhWvhWB/T5J/87P30JeOIgR84/Vpxhdrq96ja2JbxND8ubuHGMtV+K2nz1O+QetZZ8wgk4iQiBgEUrJYdfj8LhHeoASRYhfx4P4uDvQmGV+skY6Z7OtKoGmCn39kL69ezPPsmQVGOuMs1lzCqicitB6NdvDtE7OYhkZaCxuwNtyATMxAEyKsYyQlhi7QhEDXNe4e7eAn7hri7FwFU9c42JvEMi5dRO8Z6eBP355dNb679mTb0uJxORXbW+6XpQnB4f4UxbrLQDbGjx/pIxHR+c7bM3iBZKbYWI6niZk6vpQYukbd8ZoBx4Js3ETXwiw139DpTEboTUfRhCBi6jy4vwuzeQMK3XJbdug3zOH+FA/u7yQ9bXB2rsJ8xUbXgrBkgJSYzZ5Y9aaiFuJSvJGUklTMxPEluVrodizVPQz9kgi6Z6T9Ylc2mqipr/r3R27tZ7ZoU3c8Gq5PvuYgpSDZrKTccH0GslFsN8CV4AcBRrMEQjJiYDR770lC15uuCT5z/8iqGjt3j2SZyNdWrbnbBtOrxtKuVBoer1woMNIZZzxfo2r7SLnUUFqj7voEcrVbvyMRoS8TZbQrDhIeOtC1Kn5op6MEkWJX0d2sMrsSy9A40p+iM2FRc3w0TVCxPQIp2dOZ4NhQlrFcjVLDZbrQQIjwJp6KmvSmI5ydC91fri+JmBoD6Sg/driHZMRgMBvjtYt53p0tsbcrwR3DWXQtNPubusaJySJ+IDkykOL2oe1Ri2i+bK9yPYbzYZGOmvRnwpieB/d1ETV10lGTUsMlETFIRgwWKjaVhoehmcQtg7Lt4gcSXROkogb9mSh96Sj7u5McGUjyznR5WQwtcaAnuanHu1aODWUo1FwyMZMfnprH9eWydUGIsJddZ8IgCCSpqMFQR5x7R7P87ovjOM2eFFJKpAwb2A5mo5i6xp17sttmDjaSke44/ZkoJ6fDKu9CSnwZsFSeKmbpxC2DuBlaJherLjFTZ39PgsP9KWw34ORMicWaS8LSuXUgzeO3rLY8HuhJ8um7hjg+nsd2A27pS3L3nu0hRsdyVWzXByFIRswwaURKTF1juCOMy/KDAD9g2TIdMTW6EtZy5uZtgzs3XuhqKEGkUBA+MR4bynBisgiE1oyuZISPHeunNxXB0DTu39vJj84uIMuhGMjGTe4czhIzdQo1B13TiFs6P3P/CEPZOPNlm9976eJyHZqxhRozxQYfuz30xx/uT23LbtFdSSu0aFxmqVkpND9++wASGMjEEAIO96V4Z7rEC+cXiSR1qrZHvuZgaBrdyQg1xydq6dwxnOXHDvXwxOEehBAc6ivz5Ltz1Bwfy9B49GD3suhqd+4Z6WCx6nByukRfOspAJlxXM6UG08UGni8ZyITib6gjxoeO9tGVsPjTt+eYr9jLPdCSEYO9XQk+c//IFh9RexExdD55xwAvjy1i6hoBgoihL1ewjlk6g9kYnh+gVQSZmMXRgTTpmEnE1Pjg3b388NQCFdvDMsLz+2oVmPd1J9jXvf0qM+/piFFqeJSbJTESlo7tBRzoSXLXniyz5Qb5qsNMqYHrS6KmRjZuMtIVx9QF94x0cHCblgRZK0oQKRSEcR8fPtqPrgnqjo+uCdIxk0cOdtOZsLhzT4bXx4uMdiXQNI24qXHHniyGpnHbYIaPHusL4xgyUaJWeFq9djF/RVG+d2fLPHKgO2yGuk1JR03uHe3g5bH88raOuMlde7KXfk9Y/NxDoxRqDqaukYgY+FLy+kSRhhvG0kAooh7Y14mlayQiOr/4xMFV7ojD/SkO9CQoNTxSUeMKa1E7o2uCj9zWz+O3dPPodIlnT+cIpORAT5Lj4wVGO+N0JK1lN86L5xf5qfuG2d+TQNdEeCMXYrltguJKRrsSPLS/izcni8u9hT0/wNAEqahFZ8LC1DWKdY9k1CAdC8+7Rw90c6g/zcHeFIW6SyKiEzHa3w12M+i6xmA2yrszodvVMsLz8KPH+nnfLd18+61ZpJQ0vCB0+WuCx2/p4chAmrilbwu3YKtRgmgb8+iXv39FOupu8ve2mtuHM/RlIpydqxIxNW7tTy93C//AkT4O9qR4a6rIc2cXyMSs5RL+HXGTw33pK+J/qlep5ipluH07CyKAx2/pYW9XYjke61BfalWM1BLZ+KVWEveNdnJqpsJ8uUGlEVa37kxGLr1GQOQqn2HoGp2J7duSIm4Z3Dfayf7uJKdmy83WJgLjsgrLVdsjHTV59GA3qWiRfDUUk91Ji4f2d23N4LcBd+3JcmKyuOz2sQydoWyMoY4YxwYzpKIGn75rENsNmhaSS/WaNE1s67X1XtRsj9uHsvQmo0wW6kQMjQM9CToSFrf0pXj2bI5ywyO2ouDibYOZbX9tWg9KEG1jJgt1xr78ia0exo6iNxW9ZuXjka44I11xbh/O8MK5RfI1h+GOOI8c7LpqMPRoV4KxhdUZKnFLp3eHdIne0xlnT+eNtxzpz0T5C/cO8cL5RRxf4nj+qvePdl0quLgT6UxcEjYT+Trn5qur/r636Zb56LF+MjGTc/NVYpbOvaMd7FcxQ9fkgX2dzJTq/I9XpwDoTkUYzEQxdI3HD3UTt3bnba47GSEVNYAofSvczHu7Ehi6xl+8dw/Pnl1gqlCnqym6d7MYAiWIFIqbZrgjzvC91xcCdw5nmSk2ODVbRspQDH3s2MC2qTO0EYx2JRjtSvDx2wf4o+OTzJVsALqTFh840nedd+8c3n+kl1JjioVyePw9qQhPHO4BwtiYJw738sThrRzh9kEIwU/cOUQ6ZvHGeJFASixD4wNHenetGILQ+vXRY/388ZvTyyU+DvYmuWckC0AmbvLx23dHfaEbZfeuFoVig9E1wcdvH+DRA91UHI++VGRXi6GVJCMGn3twlLlyAynDwpk72Tp0OemoyV96aJS5ZjPWndZyY7MRQvD+w73cN9pBse7Sm4pe1YW72xjuiPNXH9vPbKlB3NJXubAVV6IEkUKxwWTi5q43RV+LdmvMutkoIdRaUlFzW/UD3Ax0TbRli5F2RElohUKhUCgUux4liBQKhUKhUOx6lCBSKBQKhUKx61ExRIorGMrG2PvFb12x7dkvfmCLRqRQKBQKxcaiBJHiCq4mfC4XSAqFQqFQ7CSUINoGXK0iNaiq1AqFQqFQtAoliLYB7VCRWrnRFAqFQrGTUYJIcUOs1432XlYuJaoUCoVCsdUIKeX1X6VQKBQKhUKxg1Fp9wqFQqFQKHY9ShApFAqFQqHY9ShBpFAoFAqFYtejBJFCoVAoFIpdj8oy2wKEECeAxlaPo42ISimPbfUgFAqFQrF7UYJoa2hIKe/b6kG0C0KIl7d6DAqFQqHY3SiXmUKhUCgUil2PEkQKhUKhUCh2PUoQbQ3/bqsH0Gao+VAoFArFlqIqVSsUCoVCodj1KAuRQqFQKBSKXY8SRJuEEEIXQvxzIcS0EKIihPi6EKJvq8fVLgghflMI8R+2ehwKhUKh2J2otPvN40vAzwN/GcgBvw78PvDYFo5pyxFCCOCfAn8D+K2bfLvy965GrOE9u3oOv/byOBP5+vLvf+9Dh9Yyh7DL5xHgD16d4EKuBqh5XCtV2+O3njmPH4TTsMZ53BVzOF2s87svjl+x/UNH+zg2lFm56YbnUAmiTUAIYQF/F/g7Uso/bW77GeC8EOIRKeVzWzrALUIIsZ9QBB0DLm7xcBS7kPmKvdVD2DHMl9VcrpfFqrMshhTvzbXW28I6zmnlMtsc7gJSwJNLG6SUY8AY8PhWDKhNeBg4B9wOnN/isbwnX/jCF/jCF76w1cNQtJi+VHSrh7Bj6EuruVwv3ckIurZW49ru4lrrrXcd57QSRJvDcPPn5GXbp4A9mzyWtkFK+RUp5V+VUs5s9Viux/Hjxzl+/PhWD0PRYh6/pZuIqS6DreDRg91ETX2rh7GtiVk6jx7s2uphbAv60lFuX+0aY09nnMP9qTV/pnKZbQ5xIJBSupdttwH1WKVQbBG96Sh/5ZF9nJ4rozwV66MnFeGvPLqX07OVrR7Ktube0U5GOhNcXKxu9VDang8e7ePWwTTThTqdCYt93QnCsNS1oQTR5lAHNCGEIaX0VmyPAGrVKxRbSMzSuWM4u9XD2BFETZ3bhzPXf6HiPelJRehJRbZ6GNuCoWyMoWysJZ+lBNHmsBQKP7Di3wCDXOlGU7Q5j375+0wW6ldsH8rGePaLH9iCESkUCoVivShBtDm8DpSBHwP+K4AQYi+wF3hqy0alWBOThTpjX/7EFdv3fvFbWzAahUKhULQCJYg2ASmlLYT4deBfCiEWgDnCOkQ/lFI+v7WjUygUCoVCoQTR5vFPAJPQQmQC3wb+9paOSKFQKBQKBaAE0abRDKb+B83/FJchpXxiq8egaE/qjs9CxaYjYZGM7L5LVs3xyFUcOhMWiV14/BtJw/WZL9tk4ibpqLnVw9lwltZSV9Iibqm1dDlqRhQKRdvy2sU8z5xewAskmhDct7eDRw92b/WwNo2XxxZ57mwOP5DomuCBfZ08tF/VqWkFb00VefLdeRwvQAi4czjLE4d71pW23c68NLbIj1aspYf2d/HAvs6tHlZboSqSKRSKtiRfdfjhqXm8ZoGgQEpePL/I+GJti0e2OcyXbZ4+vbDcysEPJD86m2O6eGWGo+LmqNgef/bOHI4XACAlHB8vcHZ+Z9ZQmis1eOaytfTsmQVmS40tHll7oQSRQqFoSy4s1pBXKZa41EB0p3Mhd/USZWMLu+P4N5LxxdpVe4bt1Lkdu8Y5M7agyuCtRAkihUKxYSxWHSq2d/0XNik3XPJVB+Ca8ULJ6M739Esp8YJg2YIBUHd9bM8ntQuOf6003DDe7HoNUq+1tq41txXbY7G5Lrcjlx9v3fFxvIBUM26qUHMoNS5vpLD7UGeWQqFoOYtVh2+9Oc1C2UYIuKU3xYdv68PUr/4M5ngB335rhnPzFaQMK/V+7Fg/PanIqq7W6ZjJkXX0KtoOzJUb/PEb0+SqDidnS5iahuMH1B2fqKlzuL/C4f7UNedyt/KjszleubCI60sSEZ0fv7WPAz3Jq752uCPGUEeMyfwl92MionPsst5Ynh/wnbdmOT1XRkroTlp84o5BOhPWhh5Lq7mlL8lLYxaThTqnZ8vUHJ+YqTPaFee1i3nmmufY3u44Hzs2sGt70ilBpFC0iKFs7KrFGXdjBetvn5hhoWxTrLtMF+u8MVFgqlDn84/uXXUjny/bvHh+kRfHchSqLl1Ji9mSzRsTBU7PVvj8o6P88ZszjC1UGchG+cjRvh17sX53pszr43l+eGqBmuMhBOhCMFmokYyYDGVjdCUsfnBylj97Z5ZMzGR/T4KH93dzdDC91cPfUs4vVHn+XG7596rt8ydvTvPXHt+/vF4mCzW++sI4FxdrJCMGI10x6q6PoQnuGcnywP6u5Sy+XCVcly+NLTJfsRnOxrEMjYWKw5+cmOZzD45uyXGulYuLNXxf8vSpeequT186yoHeJD94d464ZdARN5nI13jy3Tm+8foUn7xjkEcOdNOf2V2tNpUg2uYIIZ4AfgB8TEr5bSHE54HfBm6VUp7cgP3dDfyfwP2ABbwJfElK+f1W72u7cS3Rs9sqWJcbLrOlBqW6y8npEkvOi2fPLtCfifKpOwcBuLBQ5Ve+fZJC3aFU99A1wRsTeSxDp+74nJsr8850kXtHOxls9ir65pvT/OWHw5vRmbkwAPaW3hQxq/1FUrnhcna+iqEJDvYmiZo6nhfw1JkFnj49x8RinXTM4MRkgbLtETV0oqZGqeGR7DEY7ojxxkSR2VKDmuPTn4lyfqHKfNnG8QPu2pPd6kPcMs6tCIauOx6Fuoupazx7Zp5s3CITNfkX33mXYt2l7nhMlxqYmsZHj/VjaDpvTZXoTUfpiFu8NVXkd18cR9cEddfDD6BU97hjOIMmBHMlm1LD3TZp+mfmKnz9lXGeWwqibibROZ5PueFh6IKIoTNdqNPwAsYXa9Qdn8l8nZ99cISu5O7pqaYE0c7jW8DDwFirP1gIcRh4GjgO/CJgN39+TwjxaSnlN1q9T8X2w9Q1dE0wU2osiyEpJRpwZq5Mse5yIVflH/+PN5eDW11fYmpgewGaJvBl+J6FqosXSB7Y14muCSoNydOnFjgzX1mOr3n69AI/effQsmhqR84vVPnm61PLGXPPnlngE7f386/+9DQnZ0pMFur4gcQQUHV8pIS8dDF1EELj/HyN3lSUuutRscObGITZQvNlm9cu5ne1IFqyAk0X61zI1ZBSkqs6nJmrcO9oB2fnK5ybrxCPGEwU6rheQE16fPvENEPZOJomODtf5UKuykShju366JpAE4K+dASwKNRcOhMWuiawtpG78pkz8zx/Lsf4Yo2q4wNQqrvMlcJjkEg8XzbPORAC3p4qoQnBkYE0Hzrat5XD31SUINphSCnngfkN+vi/D9SAj0opKwBCiO8AbwBfBJQgUhA1dY4NpXlzsgiEga6FukvDDXjlQoF7R/N8841pSnUXP5AEMvyv4shQQDVFgwBsLwySfeFcDsvQEQJOzpS4dSCN1qwX43gBT5+e5zP3j2zNAV8HKSU/ODm3LIYAao7P//vkWS4u1ijWXaQMj6PqS1aGAzs+CAKqjsvYQpWYpRNISTJyyTrh+ZKGG7CbOTaU4eWxSyUZGq4fCm3P5/lzOfJVh1zVwag6BFJSd338AGqOTa7qsLczTs3xOL9QxXZ9hBDomiBq6syWbDIxC88P5/i2wfS2ctu+O1Oh0vDwAokQYHvhCvOCgKauxl+x6KSEsu0xka9zarakBJGiPRFhxbD/jdAqMwA8B/yny17zeVa4zIQQXwI+D/wt4J8BR4DzhBWzzwD/BniMUET9qpTy37zHEE4AZ5fEEICU0hdCnALuXv8RKnYKTxzqZbZk8523ZshVHLIxk6ip03A9fu3PznB6rkLddnGa93EBXJ4XJAEvgOlCg2LdpTtpISWct33GF2sc6ksxmI0t37TalYYbUKy7LFRschW7aXWI8s50ibmSzWLVxvEl10qMWpKHVddDCtCEYLHqYLsBqahBZ8LiYO/Vg4d3C5mYyfsP93JqtkKjad2JGDrzFRvXD/CCgKrtEcgr15ntSd6dqy55khCAEBKvKdZjpkEionPrQIrbhrLcvc0scZYhqNoeNcdfJXyAK35fIpAwX24wsbi7al4pQbS9+GfAPwJ+BXgS+Ajw72/gfb3AvwX+MaHw+ZfAV5v//k3gXwF/F/jXQogXpJQvXe1DpJT/+vJtQohO4McIxZlCAYCmCT595yAXclUKNQcBlOsOcxUHxwuoNDxW2jTeK0naDSSVhofjBmi6QCCpNS0AizWHO4Yy7OlMbPARrZ2oqbFYdTg1W6bh+rh+wNvTpeU2Cte6Ka2k4fhUGx5xy6Dh+tRdn5odZp0d6U/x+C27p3r3tTjQm2QwG6VQC0s3nM7XyFUdgiCg5l5/kuXKnzIURq4viRiS2wbS/NxDe9G07VfFuuH4zFfsG1pnK/EDyZm5Cq9dzHP3SAdz5QbzZZu+dJTuHRpXpATRNkEIkSZ0Wf26lPJ/b27+rhAiBfz167w9BvxdKeX/aH5WH/A7wK9JKX+1ue0t4CLwCHBVQXSVMRnAfwTShGJNoQDCooLfeH2Kct2l5vhUGh6lhoPthk/eN3ltxpdQ9wI0Dww9jHMo1l3SMZOFisPPtnHWTyCh7npMFxuU626YQn+TLi4fyFUdinWPpfAVTYOupMVoV3xbuXA2Aj+QfO+dWc7NVzkxVaBQ8/D9AO9mF9oKJGBpgsFMFCEEY7kq+6+Rxt9OSCmX2498960ZvvnGFN4aPKp6U/y9eH6RhYrDiaYLHODukSxPHO5tyXjbCSWItg8PE2Z1/dFl27/K9QURwLMr/j3b/Pn8im0LzZ8dNzIYIUQM+G/Ap4G/LaVUFiIFABdzVf7ZH7/D2bkKpbrLYs3BD97bCnSjBIRxNboIqDo+UVPjsVu62dMZb8Gnt56XxhbDsgLnFslXbKrriPXxgjDuQxBaL+pOwAvnc+QqNj/38F5+/EjvrsoIWsnJmRKnZsvUXD8MEA7WJ4aWiJg6thdQarhUbX/9H7gJ/D9/dpqhbAyJ5J//8UmqztrWXN0NeHu6yGLV5oF9nfSmLyUtvHaxwOH+FAOZ9k1kWAvbJ1ResdTR8fKA6ekbfH/pKttW1m2/4cuHEGIQeAr4FPCLUspfv9H3KnY25YbLV18c58REkbmyTaHm4LVIDK3El2G13TOzFd6dLvHtE9M8c3qhrartnpgs8szpsKZQoe6uSwytRBIKwwBwPUm+5vDcmQX+8PgUwXUqNO80/EDy5kSRb7w+zdm5Cq4XULX9NVlErkbN9qg6Pi+ezxOzBOfmK3znrRl+8O4cC5X2jFuTEl4fL/Cvvnuacn3t54MEak7ARL7OM2cWVlVNB1YVtdwpKAvR9mFJCPUDr6/YvqnBA0KIY8C3gSzw56SU39zM/Svam3emy5yYKlJsuNiuz0YmP0nClgq/9/IEFxZrHOxNcWKqyGfvHyET3/oaMW9Ph88gpaabbCMQAgSChYpDqe4yka8z0tWe1rKN4FtvhkIoV7GZr9jMl23cVqkhwJPh9xezNP6Pb7zDwd7kcruLExNF/vy9wwxlYyxWw/nvz0Tbwn05lqvi+VcGUd8sEnB8yUyxwfhilQO9l6rEd2yzat03ghJE24fnCFPePwN8Z8X2T2/WAIQQI8CfNX99n5Ty1c3at2J7MFWoMVdqEEha9pR+LYQI4yVA8u5MhYO9KeqOz2vj+baIbwgCydtTRRZrDt5670zXwPNls3hgONli+8X8rpm5UoOzzeKcvekIc+UGM8UGrV52DS9gYrHObNFmptTgsYPdpKImXiB5/twCMdPg3ZkyAJah8cFb+zi8xe1lhBDLcUTrRQJuANPFxrIgGuqIsa+rfRMZ1ooSRNsEKWVVCPHLwL8QQpQICzA+DvzNTRzGbxBmrP2vgCWEeGjF3xwlkBRjCzWmiw1cP2i5m+xyhAQvkNhewGLNpuGGWVfFdbgJWomUUGp4aIimcNuAfQCuH7qOOuJhe4/dwhsTRV67mGemWKdQdynVHBobFObjSxB+wELZ5ux8lTuHMwgheHe6TMy6dBt1vIA/fXtmywPd93cnODlTauk5OJmv8fLYIvfv7eTxg93bMuPueihBtI2QUv5LIUSZMNvsFwldZ38d+N2N3rcQIgF8rPnrFen3hIHa/Rs9DkX7cmKyyDszpet2Gm8VQhMYWhhgbAjBO9Ml7tqTZbijTVxGIuxjd2q2xAbpIQB0LWyG+5N3D6NpgvHFGucWqiQsnaODaeLWzrvMT+RrvHgux7szJYp1F9u7+czFm8WToYv2/HwF1/PpSFhkYtYqQQRhqv50scG+7q2xoJi6IBEx8NzWqkMfqDoez53NMVe2+fsfOkR3amcF8e+8M2WHI6X8TcLaQSv5vRV//x3ClPql378EfOmyz/g2yx1tlrc1Lt922d+rqCB8xXvw5mRx2RKyGaIoE9UxdA3Xl/RnY9heQDJqcOdw5vpv3hQk4/la2Mz2apUnW7IHSEcNDg+kMHTBj87mVjU5ffVins/c1x4xVa3k5bFFXhxbpGr7NFqRTnaDOL5kqlin1HDpiFsc6E1ie/4VIjwV3bpb60P7u/gH//04Zae1gqhm+8yV7NBVjeTXvn+av/Pjt+yomkRKECkU1+DRL3+fyUKYSTHTvMns/eK3dpVb4mZYEkG+vHbV5VYhCNs1dCej9KUjVGyfiKHxgSN9GG3QZ2qu3OC5Mzmmmz3KNtJCdD5X5yvPX+TMbAVd0xjpjLH0bFO1fV65uMgHjuyc9gvFmssfvjbJ2YUKtRbf9G8EKcNYoWTEAAlzJZv+dHR53R3sTW6pSPi9l8dpuAFBi4Opqo6P7fl0xC2QYLs+L55f5OO3D7R2R1uIEkQKxTWYLNQZ+/InAHji+V8F4Mnm74oruXUgxR8dn8DepL5a5+arRAydxapDOmZwS2+K0TbIsJJS8idvzlB1fHrT0eVK3RupEYs1l5fG8gxmo8QtfdUNOVdxNnDPm0ux7vJ/f+8UZ+bCFh2b3cFt6XsUCLxA4gaS2wbTHOpLoWmCkc44x4a21kI5X7bRNDZkbvwAig0XQxd0JiPkqjtnbYESRAqFYp0EgeTcQoXz81Vmi40Nj+WAZjqwFzBXsolaetiZXBMYbRDoWai5LFYdUhGDcsMlbhnhDWoD7962F1CsO0QMjfmyvUoQDe4Qi6bnB/zeSxd5Z7oUNsrdgpJLS209vCDA0AWpqEE6bvKx2wewjK23TALETf2KmkGtQhKe74PZGJauMZiJbsh+tgoliBQKxZoJb1LjTBbqnJgqbqo1wjQ07hjOkIlb6JrA8QLGF7e+Dk/U1NGEoDtl8fpkAdvxMTQNbwMV0ZJAlIFkfLHGYDZKJmbRnYpw7+gNFZ9ve84tVMOedl4Q1nUSG213uzqCUBSUGy49qQhPHOptGzEEsL8nwXNnF67/wjUgCFt6dCYssnGTB/d3Xfc92wkliBQKxZpwvID/+Ox5fnByDggzf9xNrJTcmbDovCxWw/G3vr1CzNLZ2x3nt56exXF9bC9g43OgAATpuMm9o51ETY1P3zXI3q7EjkmPXnLFBlJi6RqmJvCDjY9XW4mhQW86yp3DWdJRg2NDGY4OpjdvADdAImKib+BXPpSN8em7hrh1IL3c72yn0D6yVqFQbCuePbPAO1OXOsLUbH/DChAuUXjmKwD0JC3ed7Bn1d8iptY2Pc3KdQ/L0HB8SaEWpoWvZOk4WoUmwmyzQ30pMjGTiKEz0hnfMWIIQsuHEKEFrjsZoTsVIWbq106NvQlu9PtYskTWbY+YZayrNcZGYXs+lav0XVvvmhOE593ffN9+jg1ldpwYAiWIFArFGjk9VyYRCY3M5YaL6298kGvx2a+iAR+7vZ/333rJVZGJmXzqjkEixta2TfD8gFcv5Hl7uoTnBeSrzlXnpPjsV1u2T11ANmYy2pWgp1kXpitptUW2XStJRAw+defgcsHD3lSUDxzppRWHeSPfhy7CDLOq7fHqeIE3JvKbap26UaKGviFrLh01+OQdA3zsjsF1fU47o1xmCoViTVi6RkfMYLJQZ7Fq01hjV+2bZagjSqHucd/eTu7ck6Xm+KSjRstaFayVYs3la6+Ms1h1ODtf4VyuuinlB3rSEYayMUY7E0QMHUMTvO+Wnuu+dzsyWWgw0hmjZnsIAWfnK2xQm7hVGBoYmkYgJX4g8QOfi4v1DetRtx4CJKLFdR4sDQ70JPhfPnAL6ejOqmm1EiWIFArFmrhrpIP/8PRZgiBAIDYtxlXTNHqbsUOmrpGJtYcl5LmzC5QbHqau0ZW08DahfUk2ZvAX7h7iJ+8ZptzwCKTklr7UjrxpzZUbvHohT2ciwr2jBuP5GucXqhu+X02AJgSBlGhCNBvqQtzSGduE/d8sjudj6AKnhWrc1AW6pvHbz43xt544SMza+ga2G0F7XEkUCsW24649WdJRE13TMDSxaY1Fe5MRHjnYftktS0U8AQ70JOmImy2Jb7kWAsjELW4fznKwN8XdIx3cO9q5I8XQqdkyX395gnMLFaq2R8TUEUJgGdqGx7JohMI7kHI5U1DTBOmoibGR0ctrZKZot9xdmolHQMC7M2W++/ZMSz+7nVAWIoVCsWbuHe3g9GyFmuO1/LMLz3zlqnEPX/9bj/L1Fb//8i//Ml/60pdavv+bpStpUW6E86AJwaH+NIX6IpPf/89XPY4Lv/LJVb9nHv0s2cc+d8P7i1k6UVPnO2/NcKAnyS19W9thfaN47swCL5xfpFR3mSvZzJdtbh1Ik222I4kaGpWbqFh9rXV1re8jZumYusD1QytRICGmawxko9yzp/1KGswUG0z+2X9h8Zn/dsXf1rLmuhMmfiAp1V3601HOL1SXGynvNJQgUigUa+Yjtw3wZyfnOJ+rorXYZ5Z97HOrLtY6cO5XPsmzp+d55GB3y/bTKh7e381UYWK5KJ6pa9w9koX3/9wVN50Lv/JJRn/pm2vel04YSJ6I6IDg+fOLO1IQ2Z7PqxfzAKRjJh1xk3zNZapQ50h/mgM9SU7OlKg5Nx7Qf/m6gvf+PkxDoyNu0Z/V8H2JZWjNooRxHrul/dahJgSdj/8sqUd/dtX2tay5/nSEuutje2HM1lzZpi8d3TRr8GajBJFCobhpinWXJ9+d47mzOYp1h5ih0TA0NBHgNDuPtzJ+RgDRZtzCbVvcGuFa9Gei/NxDo7w9VcILAh7c18l/ePocuq6D21oLWtTS0ZpZTz3JCKU2TP9uBQ0nwF1RyuFQX4pc1UETgg/e2sfPPjDM3/ivry13vG81HTGD+0Y7SUcNelJRbM9nseqga4InDnfTm26/Ss2D2SiJiEGhvv41t1h18AKJoQm6k6EbrmJ7W57NuVGoGCKFQnFTSCn5o+OTPHVqnrGFKq4nqToBhiYwdI0lS3orHyI1wXIgZybWvjEymZjJwwe6ePyWHjqa1XzTUZNWhbkIwvRv0ZyPY4MZYpbeFj3cNoJ0zFh2jQEIIehORvjg0T6ODqb55puzxJoWnFYbLQwN7tvbyS9/6ihDHXEsQyMVDcsbDHfEOTLQXgUZl3jkQBfdqdY0l3UDiSBcd6EQhZ5UhKrdehd5O6AEkUKhuCmmig1yFYf5sg2AF8hQBAmB7fp4AS1tZqoR3vw/cKS3RZ+4Obw1VWJPZ4IH93e2pPt51NBIWDrpmEnCMnjkQDcdCYuupMX7Du3MNHshBB++rZ/4iqymwWyUB/d1Mr5Yo1R3qToehZrTMtG5REfc4lN3DjHUEeeh/V1Nl3DoknpwXycDmfbsEff4oR6yLXpo0AiLUVqGRtQ0qDk+MUsn0katSlqJcpkpFIqbImim864UPIauIQnQoyYNr9muIlifKBKAoQuSls6PH+3jJ+4cIv7Lv7yeoW8qQbMWzEAmxmhnnLrjUrbDSJfMo5+9qc/SWBKeAkvXuKU3ySfvGGCoI85gJrrlNZg2kqFsjL/62D4mC3Uihk5/s6Ho0toq1b1m2Ye1N3xd+X1YuiARMfj47QOYzSyyhw90cWwozULFoStptXUmXypq0pOMoAtYWTj+ZtccQCDDeLWQ0BX+4L7OHVf0cwkliBQKxU0xlI2RiZl0JyNMFeqYuoahazRcH09KdE0DQhfaWnqbaYQ1XkwjvPmPdie4d6STB/Z18lgbZJPdKEcH0kzm6+iaoDsVIVuyqNgNJNxUNllEB6EJgiDs39aVjDCQjXF8vMC9o507WgwtYegao12JVdv2dMRIRgwCKUlEwzYagVhbb7Ol7yNiCIazcfqzUTIxk1tXuMVSUZNUGwuhlfRmIliGRt29FGp+M2tuCU2AoQkipk46avKZ+4a5d7SzlUNtK3amzFMoFBuGpgn+3N1DPHygi4FslLil05eKkIlZRAwNKSWCppnduLmbtakL9nTG2N+TJG4Z7OlMcPtQhplSg++8tb3qnxwbyvDowW5ilt4MSo0ykIncVKxL1NTY05kgHbVIWDqD2Rh7uxIMZqNUbZ+5cmPDxt/uGLrGn7t7iP5MlIihkU1YxExtzbFEqajO4b40+3uT3L+3kw8d7Wdvd+L6b2xDHj/Yw1A2dtNtTQRh7JSpCeKmRtTUsUydQ70J/urj+/ncQ3s3Yrhtg7IQKRSKm6YzYfHT9+3hL947TKHm8NvPjuH6AW9NFZkuNpgq1AmkJBOzsFyfhusjpMS+Sm60roEhwqfQ2wbTBBIWyjZCCAayUZKR8Kn87HyFcsPdNk/pAA/s6+TukSy/8YMz7O9JIqXk6dMLvD5eeM80cQ0wDcFoZ5xs3MLQBUf6u1ZZSYT4/7f351Fy3/Wd7/98f2uv6n1Va98sy0bebYg3sE1Y7eDATBbCTSC/TJKbgQmezNyMT27mxMyZCU4g92SSE08IyVwug2NyBiZDwOAYMMZgYzC25X2RJXVra/W+1l71/fz+qFK7N0ktqbqruur1OKdPdX266vv91LdLqld/VoiHG/u/8O7mCP/lg5fx6W++TCrnk8wWODpRGuifyfv4LBzcHwl6pe7cRccJB4zb9vbwp//iCsJBb923ut22t4cnD49jrw0zOJkhnS9yun2XDehrCZMplK5WLBygLRZiY3uMvb3NNEWD/PY7dq9l9aumsf81icgFMbPSIGozwsEAV25pZ3tXnkPDs6TzPtu64jRFgmRyRaIhj2SuyOhslgNDs0xn8sTDAToSYQpFx6b2GFduaSeZLfDC8Um6miLs7GqaO5dzUDjd/+o1zHduLvyYGW/f011aZHAmTTQYIFfwyRZ9ir4jGgoQMOhpifGBK/q4Yms7s9kCu7sTPPTiEIV5/UF7NzTX9Iy7tbKnt5m/+sg1PPHGKFBar+h//vQYzx2dIJ33CRjEI0FOTmWIBI3ZbJFkeTf4prDRGo+wuSPO3e+7hEidLDboeR5/ePulfOCKjbxwfIpUrsBkMscXnugnm/dxjrn3ZCIcYEtnE54HP3tJD68OztIUDc6Nk3rrjtpbFX61KBCJyAXxrLTP1PB0lqZokE1tMa7Z3s6OcpgJBz2u2NLGprYYY7NZnjkyyUwmz9hslqMTpe0u3rGnmz29zTx/bIp80efSjS0MjKUWnKe7OUJ7Irzmr+98FH3Hs0cmODgySzQUoCkSnFvFGuCSjS3syCcIeh4ON7fQXb7o0xEPc9slvbxvX9+CY7bGwjx7dJJktsCOrgSXb25bw1dUe05OZfjpwDizmQJbO+O8fU8P4fLsp77WGD8+PMYbw7MkwkH29jXz5KExXjs5w8mpDMlcgVgwyJbOOL0tEfb0NrOhBtcUulCX9rWQyfscHp2ltyXGzRd1cWA4WX4vOsDY3B7jht2dfOjqzWxpj3NoZJaXTkzjO8elfS11ueDn6SgQich5m80W+J9PH6M5GmJ4OsvJqQyzmQIfvX47N+/pWtL10NkU4V2X9p72eFs6SuvpOOd47MAoLxybJF90bGyL8p63bFjV11JJ331liJdOTM/dLxR92hIhktkinhnv27eBQtFxeCyJYezuaeJnL+0544J3PS3r6xqspuGZDP/zp0fnWswGpzKcnMrwoas3A6X30an30ik7u5t49LURptN5RmezBANGWyxMT0uEd1+6Yd13ky3nn18a4vWhmbn7m9sTBAMBRmayBD3jmm3tfOK2i+aCJJSu087upuUOV/cUiETkvL18Ypp0rkgsFGDfplbyRZ+AZ2zqiF3QB4yZ8Y493dywq5N80V9XY2XSuSKvDM4sKAsGPLqaIvzqz/RhxlzwyeRLXTf1uC/Uanr+6NSC7kOAgbEUIzNZuk+zKOGu7iZ2diVKa+mEAvjOkVtn761zMZ3Jc2B44fswEQny3n0buGFnF+GgR1O0Pl/7+dLVEJHzls4v3FQzVJ7Wkj6HzTbPJBTw5o65XmQLxbk1iOZL5/wlwUdB6Pwsft+dkjlN+SlmpTWGADysbtfTgdK1WOZtSDpXpKNpfXQ9r7X6fTeIyKrb1b10WnLAs3U7XbkS2uKl1aMXW+5ayfnZucy1jIUD9LXW3zig89XdFFl20H2jdoethAKRiJy3ze2lHb9PregbCwd4774NNEUau/H5/Zf10VEeAG5WmhF2zbb2Kteqflza18KVW9sIlPfraI4Guf2yvrpu8TlXZsbtl/fN7QXnmfGWjS1ctaWtuhWrYY39v5aIXLDrtndw+eZWptMF2uMhfSgBXU0Rfu36bYwnc0TKs8ykcsyMWy/u4Wd2dJLMFeiIh/EqvZlZHehtifKxG7YznswRDQXmugtlebo6InLBIsEA3c0aDzOfmdFZgU1d5fRi4QCxsN53Z6L34copEInIWRV9xyuD0xyfTNORCHPZplYNCL5Ag1NpXi3PRrukr2Vu01JZXUfHU7w+NEMw4PGWjS10NVhYGBhLcmBolnDQY9+m1rmuXVEgEuHGex/h+GR6SfmmtlgValObvvH8CQ6NJOfuv3h8ig+/datC0Xl6fWiGb74wODcL6Lljk9xxeR+7expnEbxq2H90ku+9Ojx3//mjk/z8VZuWrFlUr54eGOex10fn7j9/bJJ/cc1m+lr1fx0oEIlwfDJN/723r9rxN7XF2H73g8uWP373bat23ko5MZleEIYAJlN5XjoxVdc7X6+mJ94YXTAl2jn40cExBaJVVPQdTx4aW1BW8B0/OTzeEIGoUPT58eHxBWX5Yun133nlpirVqrYoEImsstOFnuVCUi2aTOXPqVzObjK99NpN6HquqmyhuOz6WBOpXBVqs/ZS+SLZ/NIthfXv+E2aDiIiZ7SpLcZyi05vbq//v6pXy3LXbnO7ui1WUzwcXHZ9qEZ5HzdHgsuuS6ShAW9SIBKRM2qNh7hpdxfevFR08YZmLurRAm/n65aLuxdMxW+OBnnHnu4q1qgx/OwlvQvGvXUkwty4uzF2czcz3nVpL5HQmx/7XU1hrt/VGK9/JdRlJiJnde32DvZsaGZwMkN7IkRPs2ZEXYiupgi/fuN2joynANjWmZhbZFBWz8a2GL9x0w6OjKcIBYwt7fGGWr9oS0ec37hpB0fHU4QDAbZc4J6D9UaBSERWpCUaomXD0iZ3OT/BgKdtFKogHPTY3cCtm5FgQIP3T0NdZiIiItLw1EIkUiXrfTq+iEg9MTd/MQwRERGRBqQuMxEREWl4CkQiIiLS8BSIREREpOEpEImIiEjD0ywzWbfM7EUgU+161JCoc25ftSshIrIeKRDJepZxzl1b7UrUCjP7abXrICKyXqnLTERERBqeApGIiIg0PAUiWc/+ptoVqDG6HiIi50krVYuIiEjDUwuRiIiINDwFIhEREWl4CkSyrphZwMw+bWaDZjZrZl8xs95q16sWmNnnzOxvq10PEZH1SOsQyXpzD/BR4NeAMeA+4KvATVWsU1WZmQGfAn4L+LtzfLoGES5k5/m8hr+OX/xRP2OzOQD+7bv26Dqeh5GZLF96cmDu/nlex4a4hgeGZvjG84NLym/Y1cnbdnbOL1rxNVQgknXDzMLAJ4Hfdc59u1z2y8BhM7vBOfdEVStYBWa2k1II2gccqXJ1pIHlCn61q7Du5Yq6hiuVPc377XTlK6EuM1lPrgSagUdPFTjn+oF+4OZqVKgGXA8cAi4DDle5Lqvqrrvu4q677qp2NeQ0dvc0VbsK615fS5TmqNopVmJnd4Kgt7Tx56Le838fKhDJerK5fHt8UfkJYMsa16UmOOfud879hnPuZLXrstr279/P/v37q10NOY0bdnVx8YZmPDvf3jLxPOMDV2ykqylc7arUvHg4yB1XbJwLkNFQgNv29tDXGjvvYyqKynoSB3znXH5ReRaIVqE+IlIWDnq8/7I+MnuL1a7KutbTEuVXr9/OTGbxf3Oy2I6uBP+/G3cwmysQDwUIBi6sjUeBSNaTNOCZWdA5V5hXHgGSVaqTiMwTDQWqXYW60BwNVbsK64LnGS0VulYKRLKeHC3f9s37HmAjS7vRZB258d5HOD6ZXlK+qS3G43ffVoUaiUijUSCS9eQ5YAZ4B/AlADPbDmwHHqtareSCHZ9M03/v7UvKt9/9YBVqIyKNSIFI1g3nXNbM7gM+a2ajwDCldYi+75x7srq1ExGR9UyBSNabPwRClFqIQsBDwMerWiMREVn3FIhkXSkPpv535S+Zxzl3S7XrICKyXmkdIhEREWl4CkQiIiLS8BSIREREpOFpDJGISA2YzRZ47eQ0haJjT28z7Qlt31BJzjkOjSYZmsrQ1RxhV3cTgWX2wmoUzjn6x1IMTqbpaApzUU9zQ18PUCASEam6k1MZvvrMMXIFn6Lv+N6rI3zw6o1s6YiTK/i0xRWOLoRzjq8/P8jB4dm5ss3tMT509ea5EJAtFElmi7TFQniekSv4zGYLtMZCdRMUCkWf6UyBlmiQ77wyxCuDM3M/e6ppnPfs20B3UwRr0P3oFIhERKrsiYOj5Ao+J6cyHJ1IUSj6PHFwhCu2tNGRiNDdHOH9l/XRoVaj83J0PL0gDAEcm0jzxvAsF29o5omDozwzMEG+6GiOBtnYFuPwaJJcwScRCXDrxT1c1NtcpdpXxksnpnjs9VEy+SL5os9EKkdPc2kLyCPjKX5yeJyXB6fZ3dPMe97Sy+b2eJVrvPY0hkgEMLNbzMyZ2XvL9z9Wvr93lc53nZl9x8ymzGzEzP5fM+tdjXNJ7RuZyXJ0PMWTh8d4Y3iWgyOzDIyl+MGBUZ47OsHzxyb55vOD1a7mulP0HU/1j/PAU0d49eQ0U+k3N0zNF30eeXWIv/juAe5/cmDuZy+fmOb/+fbrPD0wztB0htlMgW+9eJLZbOF0p6l548kc3355iEy+tPHuRDLHoZEk6VyR0dksJybTZPIFHn9jlM8/dpCP/vef8OirQ1Wu9dpTC5HI8h4Ergf6K31gM7sUeAR4EvhloAX4L8D3zOxK51yu0ueU2lb0HS8en2IimcM5Rybv4+MoOkc8HGA2W2RwMsPFfc1csbmNcFB/y67Ed14Z4uUT02TyRSZTeSZTeS7pa6ElGuSVwWk2t8cYmspwcDTJ60OzXNLXzOtDM6RzBUZnciSzRXIFny0dcfpHk+zb1Frtl3ReDo3M4tyb9+PhAKlcgTdGZgmUe8cODM/i+46AZ8xmC/z+V1/g996V5ba9PfS0RKtT8TWmf1Uiy3DOjTjnnnTOZVbh8P8eSAN3Oue+5Zz7B+A3gUuA963C+aTG3HPPPXPf+76j4DtyRZ+i71Mo+hR8h3OQzRc5MZHi+ESKQ6OzfOflk3zxR/1MppSZz2Y6neOl41MAJMIBOhNhnHMcm0gxmcoRDhjJXJGfDozz/Nc+z8mpNI+8Okw27wNwatjQyekMzjmiofX7cRkJBvjWF/+Cgu+TzpeCUDbvc3IyzYGRGY6NzZIvOjzP8J0jV/BJZgt84YnD3P/jIzz2+ki1X8KaWL+/YZHzZCW/b2aHzCxtZt8Fti56zIIuMzO7x8z6zez9ZrbfzDJm9kr5/h4ze9jMUmY2YGafOEsV/iPwbudcal5ZtnzbGH+KNbhPfepTADx3dJL/8s1X+B8/6ufEZJpUzidTcDjAd1DwYTpbZHQ2RybvM5HKM5Mp8KODY9V9ATWs6Du++8oQ/+3RQzxxcIzvvjrEk4fGGZzKMJ0uMJHKcWwyzWtDM3zt2eOMpwqM/uDvyRSKFH3HdCZPNBSYa4Ur+o6WWIgdXU1VfmXnx/cdJ6fTfPtLf8U3njvB1/ef4OUT0wzPZDk0Osvx8QwnZ/P4vk++6JMt+BSdwzlHtlAKh08PTDA2mz3LmdY/dZlJI/pj4PeBPwEeBd4DfH4Fz+sB/hr4v4ER4LPAA+XvPwf8P8Angb80sx87555a7iDOuePAcQAziwFXA39FqXtO27s3iKPjKf7+xwM8+PwJkjkft8xjHOAc+M4xk8nz+tAM7fEwJ6dXo+FyfRlP5nj2yASz2QJbO+JcvrmNgGf8+NAYD71wkgMjMxwYngVX6iIq+qWg2ZQLcGh4lsl0YcE1zxdL44pODaSeyeRpjoa4fEsbd165iR8fGmN4Jkt3c4SrtrYRD6+Pj89njkzwdP8EAKMzWZLlFrBTTs0nK0IphZfLcvh0N0XmHjc4lebkdIY3hmeJhgJcsbmNDa319ffb+viNilSImbUAvwfc55z7g3Lxw2bWTKnb6kxiwCedc/9YPlYv8AXgL5xznymXvQQcAW4Alg1Ei/RTCloZ4Becc7NnfrisV7mCz5HxJKOzpe6uf35pkAdfOHnaMLTk+UXH80enGJvJsKUzQVdTmGu3dzCVzrO5PU5TpHH+O59I5njgJ0fIlVswDo0kOT6Z5o7LN/K/9x/j0ddHmMkUyBV9fAfJXBEPMIOpdP7U5/6yig6OjmeIhjw2t0XZ09PEN18YZDyZo+j7PD0wzo8Pj/HxW3YTCQXW5gVfgGcGJvhp/yjAkjAEnDaIF/zSYP+R2QzdTVEOjiQ5NJKce8xrJ2f40NWb6mo2WuP8CxIpuR4IA19bVP4AZw9EAI/P+/7UNIwn55WNlm/bV1ifj1L6/+dfAV8zs4865760wufKOnFiMs3/euY4zx2dZDpTms30xScGyOQLKwpDp/jA0cksRyez/OTQOO1NYX7x6s2EQ0Fu3dvN5ZvbVqP6NWf/0cm5MHTKgaFZnjw0xiOvjjCezC0JPT4s/+m/DAek8z7942k+8/Br7O1tYXdPE68NzVAsHzhffJ2P33pRTQfR774yxJefOsKR8dTZH7yMgbEUo8kc12xtJ7hoLaai73h6YEKBSGQd6yzfLh4luNI5zdPLlCXnfX8un2845x4CMLOHgSeAPwIUiOrE0wPj3HPPp/jGF/5iyc9+cPc7F9xvvfHDtN30kRUfu+BgZCbHw68M8b59fTz62gi7e5rWTVfOhZhZZgp80ff5388eZzJ15hagyR/ez9TjDywpH/iTOxbcP/X7mEqVuirT+QIB781htyPTOX50cIx3XVp7q2Xcc889c+PU5jvdazwdn9LA/p8OjBMJeuzta1nw85nM+l2KYDn1/y9HZKFTQWgD8Ny88q61qoCZ3QEknXPfO1XmnHNm9gzwG2tVD1ldL52Y4uGXhtjxro/ytj0/x1QqT7rcZTHwJ3ew7T98oyLnOTaR5thEmi0dcQanMuzqXp+Df8/F9s44B4dncc4xOptlOlMA53jl5AzOnflvkrabPrIkBJzp95ErOsaTWbKFIm3xMPFwkEQkQGssxPGJ82t5WW333HMP//I3/y0fuu9xcsXS9Tjf91yu6MgVS6EoFg6wrTMx97Pt876vB5plJo3mCSAF/NKi8jvXsA7/Fvg7M5sbsWhmYeBWYP8a1kNW0ZOHxnj++BSDUxlS2SKZZcZvVILvO4bKg6zbYqFVOUet2bexlb0bmnnpxDQvnZhmdDZLrugYnslyljx0Xgp+qVVqZCbL8EyGYrG0PtSpLVUmUzmOjqfIFoqVP/l5enZgfC4MVcJ0Os/+IxMMTqUB2NoR57odKx0ZsD6ohUgainMuaWZ/BPypmU1TmtV1M/Dba1iNTwHfBf7JzP4rEAHuAnYBP7uG9ZBV1D+axPcdWGkw7yp8TgOlv+D7x5K885JeOufNCqpnjlIImUzl5mZJjSezpDIFVid2Qtgzgp4RCXoEAsZ0Js91O9p56MXBuT3BwkGPd1/aW9VtPmazBb62/zj3PXqwosf1HUxnC0RDAX7lbVvprcPFGtVCJA3HOfdZ4HcoLYL4dUrT7lcyoLpS53+MUmtQCPh74IuUFmq8yTn3g7Wqh6yuU60H2XxxVcdaeAaeGclsfm5rhnr34vEpXh6cJhjwiIQ8hqezPHd0ktV69Y7S7DMoTeHvaY5w00VdTCTzCzZIzRV8Hp63RUY1/PDACI++OsxYsrKLdxZdeabkWKqmB5JfiPp8VSJn4Zz7HKW1g+b7h3k//wKlKfWn7t8D3LPoGA/x5jIep8oyi8tOc/4fAredU6VlXbluewfHJlI88urwkkG+rTd+uHInKk8rPzaZ4fBokuZokP7RFPFIgEs2tBAL1/7U8HM1MJ6iKRoC5zgxleHkVJoL6ZFcye+jUHRkrMjITJamSJBEJMjA2NIxRLmCz+BUhh1d1RlfMzCW4o3h2SXdZZV4zwXMGKnjBRoViEREKmg6k6dQdLx1RwcPvnCCkLc0H5/LbLKzKQKzmTxP94+TzRfZ1BYjGCg1/j/dP8EvXreF1jobW9QcCRIOeExnCwxOpSlcYD/ZSn4fPpAtOIp+kf6xFH/wjy/wi9dsXvax1WxBaYoGmcnkl5Rf6HvOA+KRQDlgr1YHcHUpEInImrrx3kc4PpleULapLVal2lROruDzrRcH5xava4uHeHZgYsEO66vFAaOzWX5yeJwNrVHetqODcDDAbLbAMwMT3Lq3Z9XrsJau3NLGDw6MMDydueAwdK6KviNf8JlK5fnhgTFu3tO9oItsZ3eC7ubqjeXau6GZVK7yXbSOUivZxtboguUH6okCkYisqeOTafrvvb3a1ai4Hx0aW7CS79eePc6JqcwZ18SpFOcckWCAgu+YTuc5NJpk74bSmjGjddbF4Zzj8FiS8dks0+m1XwfnVINfoehzbDLFh9+6hWePTjKdzrO1I85lm1rXvE7zfe/VkVV5zzlKq3yPzWbJFnyi62CV7nOlQCQiUgEHhxfuunJsMk2u4BMKGMXC6qYi5yBQ3psiHPSYTOVxzmFm9LWu/9a3+R5/Y4yn+seZyeQprsYc+7NwlBaBxIxwwKMtHubWi2unBe7JQ6NkVun9VvAdxybT3P/kAP/q5p1zG+DWi/p6NSIiVbJ48LJf/rAOLDOGaDUEAx6b22I0RUIEPcMMOpvCXLOtftaKKfqO545NApAvujVpfVvMd6WlDnzn6GwK89zRybWvxBmkV2m9KygFholknueOTnJgeOasj19vFIhERCpgcfDoa42C2ZI9tyrNA7qbI7THw1y1tZ2rt7bx/sv6+LkrNvKRt22rq1lmvnOkcgVeHZzm+FSGNcqac7zyVzTksas7wZ7eFn74xmhNLci4u2f1VioPBz1CAWMsmePIMjPs1jt1mYmIVMCe3mYCVxovHJsiV/R5+55u0rkiwzNZ8tnV+8BsT4TY29dMOBCgJRbiHXu62VflcSyrpeg7Xjs5w6GRJAXfJ2CsWStRab0nMDPCwQC9LTE6EmFyBZ+JZJ4NrbURPO+8ciMPPneCSveaBb1SIGqJllog6ylon6JAJCJSIbu6m+b2EnvgJwNcurGV7MA4s6sUiMIBwzNjIpXn1ovb+T/eto3WeH1NsT/F9x1ffuook+k8Dkeu4BMMeOT9tZlm5rvSApiRgEc0FCBSHj8TChhtNXTNXz05vSqT4p0rXYNoKMCunqaqrbO0mtRlJiKyCiLBAG2xEJ6tXhPGqS6MZKbA1Vvb6zYMARwaTTI2myXkeXQlIkSCAWKhAGs1rtcotxJ5RjTo0VJe2+n6XZ01NePqfz19vOKtZgY0RQJ0NUVojYe4dnsHWzvilT1JDVALkYhIhSWzBcaSOZ7qH2d4prJbKJwSMOhriRAMBLh4QzNv39O9KuepFalcgWS2QCpX4PhkmnzBJxIKEAp4+L6/anuYnRIodxldvqWVXV1N3La3h719LTW3p9dUuoBR2aUTjVJXITh2diV496U95fv1RYFIRKTCvvnCILOZAj3NEQ4smo5fCR4QCXl4nkd7IszeDS0UfbdmM9rW2hvDM7x8YoofHRojGvQIBzzyRZ9coUgk6FH0fXKrNEwr6EEo6GEY27sS/MI1W7hpd1fNbqQbDhqhoEe2goP5PQ+CnhENBSg6SETqsyVSXWYiIhU0nclzbKK0EreDinfpGBAMlMazpHKlQDCVzvOtFwcre6Ia8aODY3z9uUGePzZFwIzR2RyxcIDmSIi+1hiRoMdqDCPygOZIgEQkRHsszJ7eJt5xUTe9LdGaDUNQGsdW6a1DPDOaoiG6miK0xUJ1u4mwWohEZFUst0UH1Mc2HWdyqo0mVyiWZkNV+LMjHDQiwdKO6y2xEN3NEcJBjwNDs4zOZumq4Q/rc5Ur+Dw9MA6UumwSkSCxcGAuBBWdKy/QWPlzBwPG1o44Wzri9LXG5sYMPT0wwdVb22t2UcJd3U2MJXMUfZ/JCqzkHfKY2w5mQ2sUzzzqsLcMUCASkVVSr1t0nE1zNMT2rjhf+ekxCkWfSg/oMIyWaHCulSI7byG+qXS+rgJROlckX0477fEQwYBRKJYWobx0cwvjyRwnp9I0RwPMZIoVu8yhAESDHoGAcXF5C5RTcgWfdL5Ys4Ho5j1dRIIeD798kplM4YLCYoDSnnwbWk7tX2Zc1NtEJFg7g8grqTZ/oyIi69j79vVR9B2YEQsFKrqAYDzisaP7zcX3WmKlv2uDnrGxzrbpaIkFaS23zAQ8j0v7WmiNhehKRNjaEefjt+4mEgzQHA3RHKnch3Q0GMBRCp+LtcVDtERrty3htr297OhO0BINXfBaQcGAEQuFiIQ8eloiXLGllZ+9pLdCNa09tftbFRFZp6KhAJdsbGYqncez0ho66QoMcjUgk3f0NoVJ5n0CntGRCBP0jNsu6am7xfLMjHe/pZevPzdIJl8kHg5y694ePnjVJqKhADOZPL0tUfrHkviuMo1xpyJQMOCxsTVKb0uUoekMUPq9vuvS3pqfYXVRTxMb26IMTmdI5YrnNQ3fo7TtTMH5XLe9g999556K17PWKBCJiFTY8HSGoOcxk8mTyhXI+46AcUHdF0EPYqEAG1qj5H3Hf/ngZRR9x3Qmz8bWWN2FoVM2t8f5Vzfv4NhEmkjQY2NbDOcc33l5iBdPTNEcC5It+ISCHuGiT/YCLrJRCgGJSJBIMMD1uzr58Fu3cGIqQ67gs7k9RihQ+x0rl25s5drtnYyn8ownc/jncE2M0tYkTZEg3c0RmiIhCkXHyakMG1pra4mBSlMgEhGpIOccX3/uBGOzOTa0xgiYkS9miYWCNEdDDM9kyBZK207EAjCdX/44sZDRHAkxmc4TCnjEyqsjt0RDZAs+HYkwUNrHrN6FAt6ClZFfOjHNC8enALiop5nJVJ7+0SSJcIDJVJ7ZM8zBDwXAL8KpRxil1aY9K7UKNUWC7OxKcPGGZm7eU1pvZ71NBIiGAvzLqzfxjeeO0x4PMZUurGgafjRkhDyPcNCjpyU6F4rMjLFkVoFIRERWbjKV5+mBCcaSOQKesaEtVpp+HyiNJ4qFE8xm8gQDHrOZPIFCHt+92dVjQEc8xLfuuplUrsjvfOkZpsuPb42FCAY8tnXW37YJ5+LQaHLB/W2dcYq+Y1d3gh8cGCFbKJL3Ty0oWFrEMhzwyPsOM4hGvNL1dg6H0dkUpiUaorclyp7eZjoSYa7b0bGut6c4NpkmFAzQkQgzmy0S8oyC75Z0KVr5q7MpzI6uBJFQYC5cdjdH2N6ZwIy6G5+2HAUiEZEKKvqOidTCZp94OMDkvBlgnU1hEuEgR8eTTGcKBCi1UgAkwkHetrOTnpbSB9B/vONS/up7Bzj1B35bLMT/cf22NXs9tSixqHuwqynC+GyORCRIPBykuyXKbDpHuuDwfUfRlRZXpOjwXWm4dE9zmI1tcW7Z083HbtgOQDDoUfTd3Cau61l7PEwyW+DkVLa0ERlvBkQojREKBz26msL0tcW4cnMb0XApElyzrQ3ff/MavHVHB+3lFsl6pkAkIlJBoYBHd3OYoelsucSRzvts60ywq7sJz0ofVrO5At3NESZSeWazhbnBvIGAtyDw3LC7i70bWnji4CixcICbL+qu2Snfa+WKLW28Mjg9NyXfM+P9l/dxzbZ2JpJ5Xh6cKu09li2QzBUJeEY8HCQc9Ah6RsDzuGJzKxtao9yyt5vgvOtZL6t9dzSFiYUDOE4FPMBK4dx3jngoSHdLhK6mCL/y1i28ZWMb05k8W9rjtMZDHJ9MM5HMsbEtNtc9W+8UiEREKqg1HuKt2zt54cQUE8kcBmQLPls74gvG+8RCAVqiIT545SaePDzO0HSGeDjIf7xjL2/b0bngmB1NYe64YuMav5La1dUU4Zeu28rTAxPMZPJs60xw1dY2QgGPO6/ayFQmT64wSzjgEQ8XKfiOeCTIxb1NtMRC5IuOm/Z0c8Xm1rrtfpxM5Xnrjg6KvuPoeIq4QaHo8DyjtznK5vYoPsbNF3Xx7rf0LQmCm9pi627s1IVSIBIRqbA7rugjGDAGpzIEPKN9MkzPosHPO7sTOAcvn5jmxt2ddDVFeOclvWypw13EV0N3c4T37tuwpPxndnZyYHiWmXQeB7REg4QCHhOpPD0tUXZ0Jbiop4mORISNdfyB3xYPEQ8HuXF3F/2jSUZms+QLPuHyTL1EJMgNuzp516Ub6qZV7EIpEImIVFhbPMwvv3Urs9kCoYBxbCLNg88PlhZrBAq+T/9oioBnNEWCBDzjffs20NsAA1dXWyjg8bHrt9MUCdA/mpqbJr+rO8FFvU1895VhnjxU2g6kORrkX16zmbZ4/XUJRYIBbtrdxfdeG2ZndxNbO+LkCj7RcKDUUmRwfDLDbLYwt/hlozPnVmETGBFZL1btP4Dtdz94wVt3zN8P7eTf3w3ANf/6v/L43bddcP1O43z/VD7rdZxK53ljeJZQwHjh2BTDM9kFP9/cHuMXrt1ynqevOat2HVd8IOc4PJpkdDbHhpYoWzpifOGJfiYXDXi/eEMz77+sr1KnrbTzuY4LruHYbJbDo0kiIY8n3hgjtWhJgks3tvCetyxtaasjK76GaiESkZo1P/jc8uRnAOhfZsPY9aA1FuKabe0APPLq8JKfD05l1rpKdc3M2NndxM7u0v1UrrAkDAGcWKfvp5XqbIrQ2RRhOpPnOy8v876r89d/Lhp7qoKISBV0LjNrp1Fm8lRLNBggscx+Z51NjXHd46HAsquZd9bRZsAXSoFIRGSN3bi7a8FA1oBn3LS7q4o1qn+eZ9y4u4v5ywuFgx7X72yM6x4MeNywa+HsxXDQ4207O6pUo9qjLjMRkTW2s7uJj7xtK6+dnAGDvRta1EK0Bt6ysZXupgivD80SDBiX9LU01IDiyze30dMc5cDwDJFggEv6mmmONs7rPxsFIhGRKuhsinDDbnVXrLWelig9LfW9J9eZbGiN1v2eZOdLXWYiIiLS8BSIREREpOEpEImIiEjDUyASERGRhqdAJCIiIg1PgUhEREQanqbdi8gFm7/n2CmbVmkn8U1tMbbf/eCy5au4x5mI1DkFIhE5J6cLPxe6ketKnS70LBeSRERWSrvdi4iISMPTGCIRERFpeApEIiIi0vAUiERERKThKRCJiIhIw9MsM5EGZmYvAplq16OGRJ1z+6pdCRFZewpEIo0t45y7ttqVqBVm9tNq10FEqkNdZiIiItLwFIhERESk4SkQiTS2v6l2BWqMrodIg9JK1SIiItLw1EIkIiIiDU+BSKQBmVnAzD5tZoNmNmtmXzGz3mrXqxaY2efM7G+rXQ8RWVuadi/SmO4BPgr8GjAG3Ad8FbipinWqKjMz4FPAbwF/d45P19iDhew8n9fQ1/HoeIqvPH1s7v6/fdee87mODXENXzg2xXdeGVpSftXWNm65uGd+0YqvoQKRSIMxszDwSeB3nXPfLpf9MnDYzG5wzj1R1QpWgZntpBSC9gFHqlwdaVDdzRGCnlHwGyLTXJANrdFly/taY+d9THWZiTSeK4Fm4NFTBc65fqAfuLkaFaoB1wOHgMuAw1WuizSoaCjA2/d0Y+fbvtZAupsjXL2tfUHZzu4EF/U0nfcx1UIk0ng2l2+PLyo/AWxZ47rUBOfc/cD9ALbKn0Z33XUXAH/+53++queR9emKLW1s64wzMJaqdlVq3jv2dHPJhmYGpzJ0JMJs6Yhf0PEUiEQaTxzwnXP5ReVZYPl2aKmY/fv3V7sKUuPa4mHa4uFqV2Nd6GmJ0tNSmf+21GUm0njSgGdmi/8gigDJKtRHRKTqFIhEGs/R8m3fovKNLO1GExFpCOoyE2k8zwEzwDuALwGY2XZgO/BY1WrV4G689xGOT6aXlG9qi/H43bdVoUYijUWBSKTBOOeyZnYf8FkzGwWGKa1D9H3n3JPVrV3jOj6Zpv/e25eUb7/7wSrURqTxKBCJNKY/BEKUWohCwEPAx6taIxGRKlIgEmlAzrkC8O/KXzKPc+6WatdBRNaeBlWLiIhIw1MgEhERkYanQCQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeApEIlIRZnaLmTkze2/5/sfK9/euwbmvM7O8mf2fq30uEalPCkQisloeBK4H+lfzJGYWA/4HWldNRC6A/gMRkVXhnBsBRtbgVJ8B4mtwHhGpY2ohEpFzZiW/b2aHzCxtZt8Fti56zIIuMzO7x8z6zez9ZrbfzDJm9kr5/h4ze9jMUmY2YGafWGE93g38BvA7FX+RItJQFIhE5Hz8MfBp4MvAncB+4PMreF4P8NfAnwE/DxSBB4BvAt8GPgS8DPylmV13pgOZWTvw34E/Al46j9cgIjJHXWYick7MrAX4PeA+59wflIsfNrNm4DfP8vQY8Enn3D+Wj9ULfAH4C+fcZ8plLwFHgBuAp85wrPuAAeCzLGqdEhE5V2ohEpFzdT0QBr62qPyBFT7/8XnfD5Vvn5xXNlq+bT/dAczsw8AHgI865/wVnldE5LQUiETkXHWWbxcPmB5c4fOnlylLzvvenenJZrYJ+CvgD4B+MwsCgfKPvfJ9EZFzokAkIufqVBDasKi8a43O/y5KrUd/DuTLX2+Uf/ZX5fsiIudEf0mJyLl6AkgBvwT887zyO9fo/F8HFg+47gP+idJA7/+1RvUQkTqiQCQi58Q5lzSzPwL+1MymKS3AeDPw22t0/jFgbH6ZmW0vf3vEOffTtaiHiNQXdZmJyDlzzn2W0to/76PUYvMezj7DTKTmOOcYms4wPJ2pdlXWXCZf5Oh4itlsodpVqQlqIRKR8+Kc+xzwuUXF/zDv51+gNKX+1P17gHsWHeMhwBaVZRaXraAu/ef6HJGpdJ5/2n+c0dkcAN3NEe68ciPN0VCVa7b6Xjw+xaOvDZMvOjwzrt7Wxs0XdVe7WlWlFiIREWlI33t1eC4MAYzMZPn+62ux20x1TWfyfPeVUhgC8J3jp/0T9I8mz/LM+qZAJCIiDcc5x8BYCigFguGZDEfHU+w/Mlndiq2BI2MpfPfm6hZF3zE8neGbLwxydDxVxZpVl7rMRETW0LNHJskWimy/+8EF5ZvaYlWqUWMyMxKRABOpHM8dnSSdLxL0PMaSOb736jC37u2pdhVXTVMkSKHok/cdQc946cQUmbxPwXd85eljXLOtnbfvabzuMwUiEZE1lC0U+ZmdnTx67+3VrkrD27eplc/+82sMTqUBiAQDXLutneeOTXLlljbaE+Eq17DyfN9xYGiGV0/OMJstkMkX8cxoigTpbo4A8MyRCa7Y0kZrrP7HUs2nQCQiInUhV/D5af84/WMpEpEA127vWNDy5pzj2aOTHBiaIeh5TKRyJCJBYqEARd8xkcrxzy+fpC0WpjMR4lev31HFV7M6njs2yUMvneTkdJrpdIF0vjTDLBIMcHQixaa2GJdvbmU8mVMgEhERWY8efOEE/aNvjoEZGEvxvn0bGJpO89jro5yYytASC9KViDCRyvGT/nES4SCJSJCDI7Nk8kWcKwWr//a9g+zobuKm3fXTdXR0PMUXHj/M918bJp33KTqHc+A7iIaKBDwjmS0wmc7z8dsuqnZ115wCkYiIrHtjs9kFYch3jpMTaf7T11/ixRPT5ApFcgWfcNCjqylCNBRgcCqD4XDOZzb75iDjfNFnLJXj//3BYW7c1YXZmys6TKVyhIIe8XBtf3zmCj5F3weMSNB46KUhfnRwlG+/fJJUful2gZm8T8AAB4Wiz2snp7m0r5VYOLDksfWqtn+jIiIiK5Ap+EApCPWPJukfS9I/miSZLeLPe1w+55McTxP0IGyQLi7dTdj5UCh3r/1/T/Rz295eIiGP//qd1zk4kiToGdft6OATt+wmGKytydr5os8jrw7zgwMjDIylSISD7OiKMziV4SeHx5YNQ1C6BgUHs7kiA6NJ/q//+RwXb2jhmm3tvOctG9jSEV/bF1IFCkQiIrLubWiJ0hQJ8srgNIdGZxkYSzGbLZ728QUfTrc+cxEwB5GQx0Qqzz89d5zDo0mOTaTLz3X86OAYrbEQv35jbY0zeuLgGE8dHufA0CzOOcaTOV47Oc1kKsdkemUrUmeKjv7RFMPTWQbGkkxlcnz8louIhuq7tUiBSERE1r2AZ9xxRR/ff32YE5OZM4ahlfAMNrfFGJnJcGIqzeuDM7TEQphBOBggEvR4emCi5gLRgaEZxpKlxSYnU3mSuQIz6SwzWbekJexMfEqtRa+dnGE8mSOX9/mtt++isymyKvWuBQpEIiJSF/paYzRHQ2TyFxaGAIoO9h+d5JWTMxSKjmzB5+RMlkjQIx4O0BQJ0dsSrUCtKyvoGZ6VFltM5orkiv45h6H5ig5OTmf5h6eO8s8vD/Hv372HOy7fVNE614ra6vwUERG5AMlsnkzeP/sDVyDvw2y2SLbg4yjNxsoVfLJ5n0y+wIbW2gtEb9nUSiZfZHAqw3gyy8hM7rzD0HyFoiNf8Pmr7x2kWLzwwFmL1EIkIiLr2uBUmv1HJnljaJqfHBqvSACYz/HmzsG+g0DA2NqZqMnVxfNFn6ZIkL6WCCcm0xU7bsEvhczpdJ6jE2m2dzVV7Ni1QoFIRETWjVSuwMhMlvZEmJZoiANDM/zjM8d4/OAYB0dmyBQrHYdKHKUuFc+DLe1x9vQ21+RK1i+fmKY9HuLIWKrCwdAoFh2hoLdgGYJ6okAkIiLrwjNHJvjhgVGKvsMMrtrazv6jE3ztuROMTGfIVaan7LR8SlPyj02m2Nga5fLNrat7wvOQK/o8+vooo+WB1ZVS8H2m0nn2bW7lfz1znM3tMT5w5UYiwfqZeaZAJCIiNW8imeOfXzzJwFiS2WyBWCjA8ck033l5iOl0nlVqGFrAyl/ZvE+u6NM/mmJHjXUd5Qs+o7MZKp0NCz7EQsZ12zsAODaR5ieHx7n5ovpZyVuDqkVEpOYdHJnlpRNTTKTy5IuOqXSeb790cs3CEIBZafxQwDMOjyb50cHRtTnxOYgEA4QDq9Nqk8oVefS1YZwrXfCBsdRZnrG+KBCJiEjNm0znyM9LPhPJHLPZwpqFISgNqHbOlUdZG+MV7paqhILzKfqrMwus6ODoeJoXj08BjuZofXUyKRCJiEjN62kurURdKPoMTWcYns1SWOUxQ8vxXWml6ly+wMhsjsGpys3kqgTnoCW2OoO9A+Wx1FPpPNOZAteWu8/qhQKRiIjUvIt6mrl8SyueZwQ9I1ClmU4Bg56WCBta43Q3R3jw+UF8fw2bqc4iFPAIe5X/aDcgEQnSFg+xsS3GDTs7a3LZgQuhQCQiIjUvFg7wc5dvpCkSJBTwSOVWti9XpYQCEAl6REIBepsjbOuMs7E1ykymwMhsdk3rciZ9rVH6x2YqftxgwNjSEee9+zaws7uJSza2VPwc1aZAJCIi68L2zgS7epqYTOXw3ZuLJa6FQhEiASMeCtAaD7OtM46ZlQZae8Zsdm0D2ulsaYuRzFf2mAGD1miQeDhAJl9k36ZWNrfHK3uSGqBAJCIi64IZ5AuOTHnwUKV6zSZ/eP9ZHxP0IBoO0N4UZiZTYCJVSh0G/MNTR/n8Y4f44o/6qz6m6L7vH1y2fCWv8XQM8DxjOp0nEQnyzr09532sWlZfQ8RFROrMprYY2+9+8LQ/e/zu29a4RmvLOccLx6c4MDSL55VWiu5uipDKFSrWQjT1+AO03fSR0/48YBAOBYiHgtyws5Oh6SyxUIAdXQkOjcziu1JAG5vN8U/7T/AbN+0gGKhOe8Pptus422s8k0J5iFS+6Dg4kiSdL5KI1F98qL9XJCJSR84UeE4XlOrJo6+PsP/IJFAKR4dGkuzoTjCdyZfG7qzBeOZQwGgKB+ltjdISC9MSC/OBKzfy+smZJdtYpHJFjk+m2daZWP2KLaNQdHhQ8YUZp9J5wkGP6XR+Ta55NajLTEREalImX+SFY1Nz982MbR1xfN/R1xolFlz9j7CAlbqLCr5jT29pVeodXQl2diUIn+b8pytfCy2xIKtx+oLvSGaLpPMFHj84UvkT1AC1EImISE3KFnyKi6a0tyfCbO6I0xEP8rc/7OfVk5WfUTVfUyTItdva6W6JsqUjzjsv6WVnVwIz47JNrbx0YnpBHXtbovS1Vm86etDziEdCFNL5irYSNUWCtERDbOtI8NKJGfZtamNjnU27VyASEZGa1BoL0dkUZmx24YrQb9nYwjXb2ukfS3FoeOacNnWd/OH9TD3+wJLygT+5Y+G5b/wwPe/4CIlIkD0bWggHPVpjIXZ1v7l3WU9LlA9etYmn+seZTufZ2hnn+p1d5/YiK6wpGiT31D9w+MH/vuRny73GlYwraokGaY2F2NmdYHtXqSvwxGRagUhERGStvHffBr7x3CBT6dKsrp3dCa7b3kE46OGZ0RwLMZ7Mr3hYS9tNH1kSAgb+5A62/YdvLCjzDMKe4RkcGJrhLZta6EhElhxvS0ecLR21MwX93Zdu4OJ/83/xyM/+GgeGZ+dW817uNZ6JUVp7KRoKsK0zwdaOOLu6E3NjpjoSq7MadjUpEImISM3qaY7y6zduZ2QmSzjo0RYPUyj6fPFH/fz48DjxcJBgeUp4unD+o33DHnQmwngeJCJhnHOk8kXa4mFmsgUyeZ9rt7VX8JWtjut3dbL/6ASxcJBEJMh0unDOY6CDHoSDAdpiIT5x6y5yRTe3zADA5vYYO7qqM2h8NSkQiYhITTMzelqic/dfOD7F6+WxQ4lIEN852hNhegMBRmez5Is+uaJbURDwypPE2hJhNncm2NASJRTw8J0jHPTIFnxCnsete7tpXwetIls64ty4u4tUrsjJqQz5ok/qLH2KIa+0R1swYMRCQZoiATDjrdvb+ZWf2U6+6PPayRlGZrL0tETYu6Flyey6eqBAJCIi68rRiTTN0RBmpfEtAKlcga7mCJdubCGdK7D/2BS5gk+u6FMoOswgEvDwcRSKpagUDBih8npBQc+jtyXKFVvamMksXHXaDC7uXT9bVVzS18KzRybZ2BYj7/v4/umXro6HPC7qaebwWJJoKEBLNEhnU6lr8O17SgswhgIe+za1rkndq0mBSERE1pWWaJBw0GNHZ4L+sRStsRCtsRC/eO1mRmazfPP5k3jl6fIR84iFjPZ4iLZ4mEQkQNDzmMnkGZrOkswV2Xjrr7KxLcY1W9v40NVbeObIBM8dm8S50rYcN+7uXBetQ6fs6Epw9bY2XjwxRUs0RMgzut/+ETxjbsuToAdNkdKg9Wg4wO7eJiaSOVrjIQAu29zKHZf1VfeFrDEFIhERWVeu2tLOK4Mz9LREaU+Emc0WuGxTC//imi187vsHiYY8miIhcoUiRefobo5yx+V93LS7q7w5a5RsscixiRSjM1niH/kztnUk2NQWw/OMW/f2cNXWNsaTOXpaojSts1WZzYz3X7YR34dHXhvmwMlZ7H2/TipXxDlHLBygIxHmis3t3H5ZH62xEG2JEIlwgOePTbGpLcbu3uZqv4w1t75+yyIi0vBa4yF+5W1bef7YJDOZAls74lzaV+rSCnjG5ZvbGJxKM57MEQ6WusJ+5a3b5lo/SkL0NEeXPwHQFg/TFl8/rULLuf3yPnb3NvG/nznO/qOTNEeDxMJBir5PwDPecXE3N160cJmAd1xcn/uUrYQCkYiIrDutsRA3X9S9pHxjW4yZTIHdPW+2cHQ2hReFocZgZuzd0MJHfiZIJBRY8vNTIVJKFIhERCrsxnsf4fhpNtmMBJd+MEnl3Liri8GpTGnPLUrbaNzawK0eUJp5dvnmVp6ftw3KZZtaa2r9pFpgztXpLm0iImtnwX+k2+9+kP57b1/2gbfccgsAjz766AWf9HTBa1Nb7Iybwq6B852TXZEPpELRp38sRcH32d6ZILpM68g6cT7X8bTXcHgmw/B0lp7myIJlDOrciq+hApGIiIg0PO12LyIiIg1PgUhEREQangKRiIiINDwFIhEREWl4mnYvInKBzOxFIFPtetSQqHNuX7UrIXIuFIhERC5cxjl3bbUrUSvM7KfVroPIuVKXmYiIiDQ8BSIRERFpeApEIiIX7m+qXYEao+sh645WqhYREZGGpxYiERERaXgKRCIiItLwFIhERM6DmQXM7NNmNmhms2b2FTPrrXa9aoWZfc7M/rba9RBZKa1DJCJyfu4BPgr8GjAG3Ad8FbipinWqOjMz4FPAbwF/d45P16DWhew8ntMQ1/ClE1M8/NLQkvJrt7dz80Xd84tWfA0ViEREzpGZhYFPAr/rnPt2ueyXgcNmdoNz7omqVrBKzGwnpRC0DzhS5epIHdvZ1UQ4OEKu4M+VmcHFvc3nfUx1mYmInLsrgWbg0VMFzrl+oB+4uRoVqhHXA4eAy4DDVa7Litx1113cdddd1a6GnKNYOMAHrthIRyIMQFMkyLsu7aWnJXrex1QLkYjIudtcvj2+qPwEsGWN61IznHP3A/cDlHrOat/+/furXQU5T1s64vza9dvI5H0iQQ/Pu7D3nAKRiMi5iwO+cy6/qDwLnP+fqCJyTsyMWDhQkWOpy0xE5NylAc/MFv9RGQGSVaiPiFwgtRCJiJy7o+XbvnnfA2xkaTeaVMmN9z7C8cn0kvJNbTEev/u2KtRIapkCkYjIuXsOmAHeAXwJwMy2A9uBx6pWK1ng+GSa/ntvX1K+/e4Hq1AbqXUKRCIi58g5lzWz+4DPmtkoMExpHaLvO+eerG7tROR8KBCJiJyfPwRClFqIQsBDwMerWiMROW8KRCIi58E5VwD+XflLFnHO3VLtOoicC80yExERkYanQCQiIiINT4FIREREGp7GEImISEPZ1Babm3p/8tAYUJqKr/WJGpsCkYiINJT5oeeWJz8DwKP33q71iRqcusxERESk4SkQiYjUEDO7xcycmb23fP9j5ft7V+l8nyoff/HX367G+aQ6nHP0jyb5yeFxDo8mcc5Vu0o1R11mIiK17UHgeqB/lY5/FfAES9dTGl6l88kac87x9ecHOTg8O1e2q6eJn7u8DzOrYs1qiwKRiEgNc86NACOreIqrgL/XliP1a2AstSAMARwcnqV/LMWOrkSValV71GUmIlIlVvL7ZnbIzNJm9l1g66LHLOgyM7N7zKzfzN5vZvvNLGNmr5Tv7zGzh80sZWYDZvaJs5y/C9gM7F+t1yjVNzSdWbZ8+DTljUqBSESkev4Y+DTwZeBOSsHk8yt4Xg/w18CfAT8PFIEHgG8C3wY+BLwM/KWZXXeG41xdvn13OZTly+HqV8/9pUit6mmJzn1fKPqMzGQYT2bpbApXsVa1R11mIiJVYGYtwO8B9znn/qBc/LCZNQO/eZanx4BPOuf+sXysXuALwF845z5TLnsJOALcADx1muNcVb7dDPzr8ve/DnzRzBLOub8+5xcmNaVQ9NneGWdnd4JnBiZ47eQ0voO2WIjvvTpMRyJMRyJS7WrWBAUiEZHquB4IA19bVP4AZw9EAI/P+36ofDt/HNBo+bb9DMf4MqWWpG+VN6sFeKjclfafzezzzrniCuoiNeb1oRl+eGCUqXSe7uYIu7rjHBia4eh4mmzRxzN47tgU33j+BO95Sx8/e2kve3qbq13tqlIgEhGpjs7y7eIB04MrfP70MmXJed+fdV61c24AGFjmR18HbqPUcrTcz6WGjcxk+dYLJ/HLU+uPTaT4xnMnGJrJkCv6pHNFMvkingczmTxT6aO8cHySf33rbi7b1Hba474+NMMbw7NEQx6XbWqju7m+WpY0hkhEpDpOBaENi8q71qoCZna7mf3CMj+KlW/H1qouUjmvD83MhSGA8WSOfNEnlS2SzhVI5ooUHeSLkM77nJzO8PTABP/5G6/wdP/Essd8/I1RHnx+kNdOzvDc0Sm+/JMjDE6l1+olrQkFIhGR6ngCSAG/tKj8zjWswy8C/5+Z9ZwqMLNAufwnzrnZ0z5TatbipYUMwzOj6ByzOX/J43N5n4lknlcGp/n9rzzHN59f2EiZLRR59sjCoFTwHU+dJjytV+oyExGpAudc0sz+CPhTM5umtADjzcBvr2E17gX+BaVxQ/8ZKAD/BriUUpeZrEOX9rXwdP8EBb/UStTZFObEhDGTKSx5rCt/+b7D5QoMzzjufegVDgzPcNXWdm7c3UU6VyRfXNoDO53Or/IrWVtqIRIRqRLn3GeB3wHeR2nczntY2YDqSp3/FeAmSuOWPkdpQHcQeKdz7vEzPVdql+9gc3uM2UyBXLHInt5mPnDlpgXdaMsp+DCbLTIyneHbL5/kS0/2c/+PB2iOBGmNhZY8fmtHfLVeQlWohUhEpIqcc5+jFEbm+4d5P/8CpSn1p+7fA9yz6BgPAbaoLLO47DTn3w/cfi51ltp1dDzFPz57nKLvaIoGMYOLepr43qvD+P7Z9y9zQLrgODaRZmQ2x9HxNKGA8e639PKN5wdJ50qTDje2RXnrjo5VfjVrS4FIRESkTvzo4BjFecGnUHT82bdf49mBcZbp9TqtyXSBQKbARDLH/U8ewTn42A3bOTmVIRoKsKE1evaDrDMKRCIiInViPJWb+z5XKPLo6yMcGJqhsHQs9VkVHfhFx5HxJI+9PsLevhau3nqmZa3WN40hEhERqRMb20orJjjn+M4rQ7xxnmHoFAckM0VePjHNDw6s5h7D1adAJCIiUidu3t1FPBzg6YEJBsZS5MthaPKH95/3MYvAaDLH914d4bWTM5WpaA1SIBIREakT7Ykw1+/qZDKVIzdv0NDU4w9c0HF93+H7Po+9PrKiwdnrkQKRiIhIHXmqf5wj48mzP/CcOGazBV48PsXITKbCx64NCkQiIiJ1ZP/ABDPZCxg4tAzfd+SLjsGpNF9/fpB8sbLHrwUKRCIiInXkjZHK77hSdDCZyjOZzrP/6CQHhupvVxdNuxcREakThaLPC1//O44/8j+W/GzgT+5YcL/1xg/TdtNHVnRcAzwP0rkiw9MZBsaTXLqxpRJVrhkKRCIiInXiiYNjbH7nrxK8buGewQN/cgfb/sM3zvu4RQf5os9spoDfAtFg4EKrWnMUiEREROrE/U8OkCpvr1FJDsgXHIbPjq4427sSFT9HtWkMkYiISJ14emCC2UzlA1Gpy8xoigS5Zls72zvra2NXUCASERGpGwXfZzXmfzkgYLCrJ8Htl23E7Kz7Bq87CkQiIiJ1Ih5efiRM640fvuBj+87hO2iOhS74WLVIgUhERKRO7OyOL/vBvtLZZGdyahzR1LwNZOuJApGIiEid2NXdTDhY+e4szyAc9AitwrFrhQKRiIhInbh+Vyc9zdGKH9cz6ExEuLSvhdZ4uOLHrwUKRCIiInWiuzlCW6KyY3wCBm3xEO/Y080n33lRRY9dS7QOkYiICLCpLcb2ux9ctvzxu2+rQo3O3fPHptjVleDFY9MVmW0WDXrcccVGNrXFMIPjk5m6bSFSIBIREYHThp7lQlKtyuSLvD48W1o4yF348ToTIWKhUmeSc/DDN0bYu6EZz6u/sUTqMhMREakTuYLP8HQWV4Ew1NscZnNHgjeGZ+d2t09mi6TzlV/4sRYoEImIiNSJSNAjFPC40HUTuxIh9m1qA8B3MJXOA9ASCxEP198+ZqBAJCIiUjcSkSC7uuNcSI+WAZ1NEcIBj3CgFBOCnhH0jFsv7q7LVapBY4hERETqxlVb23nmyCTRoEcy55/zMKJI0DAMs9LeZVdsacMB77qkh929zTRF6jc21O8rExERaTD7NrVy7fZ2jowleenENA5H0Ct1oWXyRQpnmHoWCZQ2b42FAuzoTLCnt4krt7bz1u0dxOq0m2w+BSIREZE6csXmNk5MpjGDkZkcBb+UglqiQZLZAsnc0mAU8iDgecRCQW6+qIufu2IjN+zuqkLtq0eBSEREpI7s3dDMSyemGJnJ4lmSou/oSIQZT+YZS2Yp+FmsUKToQyBgtESDbG6L0dsa5ept7fz6jTuIBOu/RWgxBSIREVnXbrz3EY5PppeUb2qLVaE21RcMeFy5pZ2BsRSzmQJN0RBvv6iTV07OMDqTZf/RSUZmsqRzRYIBoy0eJhENsbO7CbCGDEOgQCQiIuvc8ck0/ffeXu1q1IzR2SwPvXiSUMDjot5mAI5MpOlsihD0PN4aDHBkLMmB4VlCQWNnV4KtnQkiwQDbOxNVrn31KBCJiIjUkddOzuAvWpkxm/e5fHMrr52cBaAjEeZnL+0lX/SZShcA6GuNcvNFjTVuaD4FIhERkTpyumWCupoi3HBDF0MzGSLBAB2JMM45hmeyeGZ0N0fWtqI1RgszioiI1JFLNrQQWLQyYzwcYGdXE55n9LXG6EiUNmg1M3pbog0fhkCBSEREpK60J8LceeVGeluihALG1o44H7x6E+GgPvLPRF1mIiIidWZbZ4JtDTxA+nwoLoqIiEjDUwuRiIjIGWxqi7H97geXLX/87tuqUCNZDQpEIiIiZ3C60LNcSJL1S4FIRETWBa1ILavJ3KLFm0REREQajQZVi4iISMNTIBIREZGGp0AkIiIiDU+BSERERBqeZpmJiEhNMLMXgUy161FDos65fdWuRKNQIBIRkVqRcc5dW+1K1Aoz+2m169BI1GUmIiIiDU+BSERERBqeApGIiNSKv6l2BWqMrsca0krVIiIi0vDUQiQiIiINT4FIREREGp4CkYiIVI2ZBczs02Y2aGazZvYVM+utdr1qgZl9zsz+ttr1aBQKRCIiUk33AB8Ffg14O7AZ+Go1K1RtVvKfgN86j6c7fS34WjEtzCgiIlVhZmHgk8DvOue+XS77ZeCwmd3gnHuiqhWsAjPbCfwdsA84UuXqNBS1EImISLVcCTQDj54qcM71A/3AzdWoUA24HjgEXAYcrsQB77rrLu66665KHKquqYVIRESqZXP59vii8hPAljWuS01wzt0P3A9gZhU55v79+ytynHqnFiIREamWOOA75/KLyrNAtAr1kQamQCQiItWSBjwzW9xbEQGSVaiPNDB1mYmISLUcLd/2zfseYCNLu9HkAjx7ZJLtdz+47M82tcV4/O7b1rhGtUeBSEREquU5YAZ4B/AlADPbDmwHHqtarepQtlCk/97bl/3Z6YJSo1EgEhGRqnDOZc3sPuCzZjYKDAP3Ad93zj1Z3dpJo1EgEhGRavpDIESphSgEPAR8vKo1koakQCQiIlXjnCsA/678JfM4526pdh0aiWaZiYiISMNTIBIREZGGpy4zERGRBrapLbbsTLNGm46vQCQiItLAThd6Gm06vrrMREREpOEpEImIiEjDUyASEZG6Y2a3mJkzs/eW73+sfH/vKp2v2cz+3MyOm1nKzPab2S+sxrlkdWgMkYiINIIHgeuB/kof2Mw84FvAxcDdwBHgo8A/mFnKOddYg3HWKQUiERGpe865EWBklQ7/YeAG4IZTW46Y2XeAncAdlMKY1Dh1mYmIyLpmJb9vZofMLG1m3wW2LnrMgi4zM7vHzPrN7P3l7q2Mmb1Svr/HzB4ud30NmNknzlKFXwB+PH//NVdyg3Pudyr+gmVVKBCJiMh698fAp4EvA3cC+4HPr+B5PcBfA38G/DxQBB4Avgl8G/gQ8DLwl2Z23RmOcyXwvJn9tpkdMLOCmb1sZh88r1cjVaEuMxERWbfMrAX4PeA+59wflIsfNrNm4DfP8vQY8Enn3D+Wj9ULfAH4C+fcZ8plL1EaE3QD8NRpjtNDqWvsNuA/AuPAJ4Cvmtmdzrmvn+fLkzWkFiIREVnPrgfCwNcWlT+wwuc/Pu/7ofLtk/PKRsu37Wc4RhjoBt7nnPuyc+5hSq1LrwL3rLAeUmUKRCIisp51lm8XD5geXOHzp5cpS8773q3gGDPAq865N+ae5FwBeBi4sjwLTWqcfkkiIrKenQpCGxaVd61hHV4HIsuUh4CMc85fw7rIeVIgEhGR9ewJIAX80qLyO9ewDt8E9pjZNacKzCwMvBd4bA3rIRdAg6pFRGTdcs4lzeyPgD81s2lKa/7cDPz2Glbjv1JaiPFrZvZ/A8PAJ4EtwC+vYT3kAqiFSERE1jXn3GeB3wHeB3wdeA9nn2FWyfNPAjcC3wE+C3wVSADvdM6dbmaa1BhzbiXjxURERGQdWPKhfsstt/DkoTEyR144pwNtv/tB+u+9vWIVqxJb6QPVQiQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeApEIiIi0vAUiERERKThKRCJiIhIw1MgEhERkYanQCQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeApEIiIi0vAUiERERKThKRCJiIhIw1MgEhERkYanQCQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeApEIiIi0vAUiERERKThKRCJiIhIw1MgEhERkYanQCQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeMFqV0BEREQq48Z7H+H4ZHpB2clDY0SCgSrVaP1QIBIREakTxyfT9N97+4KyW578TJVqs76oy0xEREQangKRiIiINDwFIhEREWl4GkMkIiIiS2xqi7H97geXLX/87tuqUKPVpUAkIiIiS5wu9CwXkuqBusxERESk4SkQiYiISMNTIBIREZGGZ865atdBREREpKrUQiQiIiINT4FIREREGp4CkYiIiDQ8BSIRERFpeFqYUUREpE6Y2YtAptr1qCFR59y+lTxQgUhERKR+ZJxz11a7ErXCzH660seqy0xEREQangKRiIiINDwFIhERkfrxN9WuQI1Z8fXQStUiIiLS8NRCJCIiIg1PgUhERGSdM7OAmX3azAbNbNbMvmJmvdWuVy0ws8+Z2d+e7XEKRCIiIuvfPcBHgV8D3g5sBr5azQpVm5X8J+C3VvJ4rUMkIiKyjplZGPgk8LvOuW+Xy34ZOGxmNzjnnqhqBavAzHYCfwfsA46s5DlqIRIREVnfrgSagUdPFTjn+oF+4OZqVKgGXA8cAi4DDq/kCWohEhERWd82l2+PLyo/AWxZ47rUBOfc/cD9AGa2oueohUhERGR9iwO+cy6/qDwLRKtQn3VJgUhERGR9SwOemS3u9YkAySrUZ11SIBIREVnfjpZv+xaVb2RpN5qchgKRiIjI+vYcMAO841SBmW0HtgOPVadK648GVYuIiKxjzrmsmd0HfNbMRoFh4D7g+865J6tbu/VDgUhERGT9+0MgBHypfPsQ8PGq1mid0eauIiIi0vA0hkhEREQangKRiIiINDwFIhEREWl4CkQiIiLS8BSIREREpOEpEImIiEjDUyASERGRJczsFjNzZvbe8v2Ple/vXYVz9ZePvexXpc+3HC3MKCIiIivxIHA90L8Kx/4gpc1o53sL8Hngv63C+ZZQIBIREZGzcs6NACOrdOxn5983syjwd8BTwF2rcc7F1GUmIiLS4Kzk983skJmlzey7wNZFj1nQZWZm95S7ut5vZvvNLGNmr5Tv7zGzh80sZWYDZvaJc6zSvwf2AL/lnMtX5lWemQKRiIiI/DHwaeDLwJ3AfkrdVWfTA/w18GfAzwNF4AHgm8C3gQ8BLwN/aWbXraQiZtYL3A18zjn33Lm8iAuhLjMREZEGZmYtwO8B9znn/qBc/LCZNQO/eZanx4BPOuf+sXysXuALwF845z5TLnsJOALcQKkL7Gz+DaUNau89x5dyQdRCJCIi0tiuB8LA1xaVP7DC5z8+7/uh8u2T88pGy7ftZzuQmYWA/xP4knPu2ArPXxEKRCIiIo2ts3y7eMD04AqfP71MWXLe9+cybf6d5fp88RyeUxEKRCIiIo3tVBDasKi8a60rAnwAOAE8ttYnViASERFpbE8AKeCXFpXfWYW6XA884Zxbk8UY59OgahERkQbmnEua2R8Bf2pm05QWYLwZ+O21rIeZecA+4J/W8rynqIVIRESkwTnnPgv8DvA+4OvAezj7DLNK66TUUDOxxucFwKrQKiUiIiJSU9RCJCIiIg1PgUhEREQangKRiIiINDwFIhEREWl4CkQiIiLS8BSIREREpOEpEImIiEjDUyASERGRhqdAJCIiIg3v/w8CtY/+Gj9MAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 49 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we visualize the posterior predictive distribution using a pairplot:\n",
    "pairplot(x_pred.squeeze(), \n",
    "         upper=\"scatter\", \n",
    "         points=xo,\n",
    "         points_offdiag=dict(marker=\"+\", markersize=10),\n",
    "         hist_diag=dict(bins=10),\n",
    "         points_colors=[\"k\"],\n",
    "         limits=[[0, 1]]*7,\n",
    "        );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Result:\n",
    "The posterior predictive samples contain the observed data; that's a good sign. We continue with the calibration analysis using SBC. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation-based calibration (SBC)\n",
    "\n",
    "To perform SBC we repeat the inference many times using different (simulated) data points and check where the underlying parameters are located in the inferred posterior, i.e., we calculate their rank under the posterior.\n",
    "\n",
    "If the variances of the posterior are well-calibrated then these ranks must be uniformly distributed. I.e., if the ranks are not uniform there is something off! \n",
    "\n",
    "We can use the SBI package to perform SBC. \n",
    "To run SBC we need simulated data. Above we loaded this data from disk. In practice, one has to simulated extra data with parameters sampled from the prior. It must not be the same data as used during training but unseen data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fceea77e0e9c429c88b6cd9583acecf2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Running 1000 simulations in 100 batches.', style=Progress…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "54deef1012f74334bfe8b0fb4dc8bcec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Running 1000 sbc samples.', max=1000.0, style=ProgressSty…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "num_sbc_points = 1000\n",
    "num_posterior_samples = 1000\n",
    "\n",
    "# we simulate with parameters sampled from the prior.\n",
    "thos, xos = simulate_for_sbi(batch_simulator, prior, num_simulations=num_sbc_points, \n",
    "                            num_workers=num_workers, \n",
    "                            simulation_batch_size=10)\n",
    "\n",
    "\n",
    "# run SBC using sbi package\n",
    "ranks, dap_samples = run_sbc(thos[:num_sbc_points], \n",
    "                             xos[:num_sbc_points], \n",
    "                             posterior, \n",
    "                             num_posterior_samples=num_posterior_samples, \n",
    "                             show_progress_bar=True\n",
    "                            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq8AAAFGCAYAAABEwnZKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeIUlEQVR4nO3deZglZX0v8O8PBtARWcIqBhGVy76oCEHDIhqDS9SIIpcrigE0ChgSr0GCy6gMIkpcwhqW5BLmQoyIQIKJYthEvAJRY9hiXEAlgoA4IouOU/ePOoNNz0xPz/T6dn8+z1PPTNepOlX1nvp1f8973qpTXdcFAABasNpU7wAAAIyW8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRjpcLrfvvt1yUxmUy/maYt9WoyLTVNW+rVZFpqWq6VCq/33nvvyiwOTCH1Cu1QrzB6hg0AANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRjzmRu7Oabb85DDz00mZuE5Zo7d2623377qd6NaeuaV71qhcvsdcklk7AnwHhYUU2rZ1oxqeH1oYceyjrrrDOZm2QG67ouVbXK6y9cuHAc9wYAmAyTGl5nszve8Y5VWm+LT31qnPdkfPzwhz/Mi1/84sybNy8HHnjgpG//a1/7Ws4666ycddZZo15nwYIFOf300/PlL395AvcMAJhIxrzSpAsvvDDf+c53Rr38ZZddlg9/+MMTuEcAwGTQ88qMdvfdd+cTn/hELr744qy77rpTvTsAwBjpeZ2h9t1338yfPz+HH354nvOc5+Soo45Kktx111059thjs/fee2fHHXfM7rvvniOPPDI/+MEPHlv33e9+dw4++OBcdtllefnLX54dd9wx++23Xy6++OIRt3nyySdnu+22y0UXXbTcZW699db80R/9UXbbbbc8+9nPzoEHHpgrr7zyccvcfffdOeaYY7LHHntk5513zgEHHPC4j/oPPvjgXH755bnrrruyzTbb5LOf/exyt/fxj388N910U0477bS88IUvHHH/AYDpT3idwS644II85SlPyamnnpqDDjoojz76aN74xjfmlltuybHHHpuzzz47RxxxRK6//vr8xV/8xePWve2223LKKafkLW95S84444xsttlmOfbYY3P77bcvc1uf+tSncs455+T444/P/vvvv8xlHnzwwRx66KFZvHhxPvaxj+WUU07JuuuumyOPPDLf/e53kyQPPPBADjrooNxwww358z//83zyk5/Mpptumre+9a25+uqrkyTvf//78/znPz8bbbRRLrzwwuyzzz7LbYNDDz00n//857PvvvuuQgsCANONYQMz2HrrrZf3ve99mTOnf5lvv/32bLrpppk3b16e9axnJUl233333HnnnVmwYEEWL16c1Vbr38/8/Oc/z/nnn5+tt946SbLllltm3333zZVXXvnYvCXOPPPMnHHGGZk/f37+8A//cLn7853vfCf3339/9t9//+y1115Jkp122imnnHJKfvWrXyVJzjvvvNxzzz259NJLs+WWWyZJ9tlnnxxyyCE58cQTs/fee+dZz3pW1ltvvayxxhrZZZddRmyDrbbaaiVbDQCYzoTXGewZz3jGY8E1Sbbeeuucf/75Wbx4ce68887ceeed+e53v5t/+7d/S9d1WbRoUdZcc80kyZOe9KTHhdRNN900SZa6T+9FF12Ub33rW9lrr71GDK5JHyQ33HDDvPe97821116bPffcM3vuuWeOPfbYx5a57rrr8sxnPjObb755Fi1a9Nj8fffdNyeccEJ+9KMf5alPfeqqNwoA0DThdQbbcMMNl5p33nnn5cwzz8x9992X9ddfP9tuu22e+MQnJunvm7rEknlLLOmRHbpM0n/xxN57752rr7461113XV7wghcsd3/mzp2bBQsW5LTTTstVV12VSy+9NHPmzMmLXvSizJs3L+uvv34eeOCB3HHHHdlhhx2W+Rx333238AoAs5jwOotcfvnlOeGEE3LUUUflgAMOyEYbbZQk+ehHP5qbbrpplZ7zbW97W972trflta99bd773vfmsssuy5Oe9KTlLr/FFlvkIx/5SBYvXpxbbrkl//Iv/5Jzzz03a6+9dubPn58nP/nJ2XnnnfOe97xnmesvGUoAAMxOwusscsMNN2SttdbK29/+9se+mWrRokW57rrrkiSLFy9e6efcaKONMmfOnHzwgx/MgQcemJNPPjnve9/7lrns9ddfn3e+8505++yzs91222WHHXbIDjvskGuvvTZ33XVXkmS33XbLZz/72Wy++eZZb731Hlv3nHPOyde//vWcdNJJSX7TEwxAb0Vfwf6jDTYYcf0bbrhhvHeJWWwiv4JdeJ1Fdt5551xwwQU5/vjj85KXvCQ//elPc/755z92B4GHH354qeECo7XTTjvloIMOyoIFC/Kyl70su+6661LL7LjjjlljjTXyrne9K0cccUQ22GCDXH/99bntttvyoQ99KEny5je/OZdddlne9KY35bDDDstGG22U6667Luecc05e85rXZO7cuUmSddZZJ/fdd1+uueaabLPNNtl4441XsVUAZoYVfQX7WkOuI1gWX9/OeJrIr2DXfTWLvPrVr8473vGOXHnllTn88MNz0kkn5WlPe1pOPfXUJP1Xro7F0UcfnU022STHHXdcHnnkkaUeX3vttXPuuedmyy23zPz583PYYYfliiuuyPvf//687nWvS9L35F5wwQXZaqutMn/+/Lz1rW/NFVdckaOPPjof+MAHHnuuAw44IJtttlmOOOKIfO5znxvTfgMA7ajhF+CMZNddd+1uvPHGVd7YDTfc4J0d08bChQvzvOc9b6xPU+OxLxNhrPV6zatetcJl9rrkklV+fpgCM7ZekxX/jb15/vwR19/+uOPGtH0Yahz+xi63XvW8AgDQDOEVAIBmCK8AADTD3QaAZbpnvfXyq9VXH3EZt9ZhMk3krXeAdgivLJOB/fxq9dXdWodpZSJvvQO0w7ABAACaIbwCANAM4RUAgGYY89og41EBgNlKzysAAM0QXgEAaIZhAwDAmK1oSFtiWBvjQ88rAADNEF4BAGiG8AoAQDOEVwAAmuGCrZU0GQPSR7MNAIDZSM8rAADNEF4BAGiG8AoAQDOEVwAAmuGCLQDGbEUXmvpmJWC8CK8AwAq5Ew7ThWEDAAA0Q3gFAKAZhg0As5qxmgBt0fMKAEAzhFcAAJph2ABTYjRXrfq4FgAYTs8rAADN0PMKADCLtH6hqp5XAACaIbwCANAM4RUAgGYY8woAMIOM5o4+LdPzCgBAM2Zcz2vrV9AxflZ0Lmx+1FGTtCcAwHjR8woAQDOEVwAAmiG8AgDQjBk35hVgMhlnDzC59LwCANAM4RUAgGYYNgBMGB+pAzDehNcJ4A82AMDEEF4BABox07/6dTSMeQUAoBnCKwAAzZjUYQP/dfrpWWvRouU+Ph3GguqOBwCYvvS8AgDQDOEVAIBmCK8AADSjuVtlGZM6PsbajtPhXrbOBQCGmg5/m2aC6d6Oel4BAGiG8AoAQDOEVwAAmtHcmFcAgFUxmmslpno8Jyum5xUAgGboeQWmzGTcMWKqe1Gm+1W7AK3R8woAQDP0vAIAk/JJyFg/iXB/bxLhFWBE/lhqA2B6MWwAAIBmCK8AADTDsAEAmAH+6/TTs9aiRVO9GzDhhFeAKTQdbprewoU6AEsYNgAAQDP0vM5ArgwGACbKVH+SoucVAIBmCK8AADRDeAUAoBnGvDIhjLsFYGXNhL8dUz0edDYQXoFVNhP+0LTAH0OA3zBsAACAZgivAAA0Y1oNG5gO3/ICAMD0pecVAIBmCK8AADRDeAUAoBnCKwAAzZhWF2wBsPLcBxaYTfS8AgDQDD2vAACj5JabU094nQJOfACAVWPYAAAAzRBeAQBohmEDAADThKGFK6bnFQCAZgivAAA0Q3gFAKAZwisAAM1wwRbTlkHrAMBwel4BAGiG8AoAQDMMGwCY4QzBAWYSPa8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzhFcAAJrhVlkAAANuLTf96XkFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADN8PWwwIzmqx6B6cTvpLHT8woAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzJvUbttb49a/z6Bxf6sX0MHfu3KneBQBgJU1qktz4gQcmc3Mwou23336qdwEAWEmGDQAA0AzhFQCAZgivAAA0w9VTADADuCia6WLhwoUTelG0sxwAZgAXRTNdPO95z5vQ5zdsAACAZgivAAA0w7ABYJmMn2MyLVy4cIXL+GIRIBFegeUwfo7JNNFj5ICZw7ABAACaoecVgCl3zateNeLje11yySTtCTDd6XkFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzhFcAAJohvAIA0AzhFQCAZgivAAA0Q3gFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzhFcAAJohvAIA0AzhFQCAZgivAAA0Q3gFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzhFcAAJohvAIA0AzhFQCAZgivAAA0Q3gFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmiG8AgDQDOEVAIBmCK8AADRDeAUAoBnCKwAAzRBeAQBohvAKAEAzhFcAAJohvAIA0AzhFQCAZgivAAA0Q3gFAKAZwisAAM0QXgEAaIbwCgBAM4RXAACaIbwCANAM4RUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDOEVwAAmlFd141+4aqfJLljDNvbMMm9Y1ifnnYcH+PRjvd2XbffeOzMeFOv04q2HB9jbceZXK+J82y8aMfxMWH1ulLhdayq6sau63adtA3OUNpxfGjHkWmf8aMtx4d2HJn2GR/acXxMZDsaNgAAQDOEVwAAmjHZ4fWvJ3l7M5V2HB/acWTaZ/xoy/GhHUemfcaHdhwfE9aOkzrmFQAAxsKwAQAAmjEp4bWqVq+qD1fVf1fVg1X1maraZDK23Yqq2r6qumVMvzt4/CVV9Y2qeriq/r2qXjps/Y2r6tNV9UBV3VNVH6mqOVNzNFOjqs6sqrOHzRtzu1XVn1bVHVX1UFV9saq2mozjmSrqdXTU7Nio1/GhXkdHvY7dtKrZrusmfEryoSR3Jfm9JM9J8tUkX56MbbcyJXl9kp8k2XTYtEaS7ZI8kuS4JNsM2vPRJNsPWf/aJNck2TnJy5Lck2T+VB/XJLVdJflgki7J2UPmj7ndkhyaZGGS1ybZMcmlSb6TZK2pPu4JbE/1Orp2UrOr1m7qdXzbU72Orp3U66q33bSr2ck46DUHO3bIkHlPHzTC86f6RZku0+BFv3o5j52Z5Kph865M8teD/+8xaM8thzz+pkG7z9hf2oPjfMagLZbc4HtoYY253ZLcnmTekMfXTvLzJAdN9bFPUHuq19G3lZpd+TZTr+Pbnup19G2lXlet3aZlzU7GsIFdkjw5yVVLZnRd9/0k30+y5yRsvxU7JLl1OY/tmSHtN3BVftN+eya5o+u67w17/Mnp238m2yPJd9O/Y/vesMfG1G5VtXGS/5HHn7sPJrkxM/fc3SXqdbTU7MpTr+Nrl6jX0VKvq2Za1uxkhNffHvz7o2Hz70qy+SRsvxU7JNmiqr5aVT+uqiuqarfBY7+dkdtveY8nM7yNu65b0HXdoV3X/XgZD4+13WbjuTsbj3lVqdmVpF7H3Ww85lWlXlfBdK3ZyQivc5Ms7rruV8PmP5rkCZOw/Wmvqp6Yvmt+3STvSvLK9C/e1VW1bfo2fGTYakPbb6nHB+3dZXa38Vjbbe5g9kjPMdOo11FQsxNCva489ToK6nXCTFnNTsaVcg8nWa2q5nRdt2jI/LWS/GIStj/tdV33cFWtn+TRruseTZKqOiTJc5O8PX0brjVstaHtt9TjVbVG+kHWs7mNx9puDw9ZZ3nPMdOo11FQsxNCva489ToK6nXCTFnNTkbP6w8G/z5l2PzNsnRX8azVdd3CJUU1+HlxkpvTd53/ICO33/IeT2Z3G4+13WbjuTsbj3mVqNlxp15X3mw85lWiXifElNXsZITXb6a/cmzvJTOq6unpr4i8ZhK2P+1V1XOramFVPWfIvNXTDwS/OcmXM6T9Bl6Y37Tfl5M8o6o2H/b4z5N8Y4J2uwVjareu6+5J8u08/txdO8mumbnnrnodBTU7IdTrylOvo6BeJ8zU1ewk3WrhxCQ/TrJffnMfuqsmY9stTOmHb/x7+ivsdk+yfZLzktyfZOP0V/n9MskH0t9L7YPpu9u3HaxfSa5P8pVB+740yd0ZcvuJ2TClv2Jx6G08xtxuSf44yYNJDkw/4P/SJP+ZZM2pPt4JbEf1uuI2UrNjb0P1Oj7tqF5X3EbqdXzacdrU7GSeOCcnuTfJz5L8fZINp/qFmE5TkqcmWZD+Br6/SPKFJDsMefzl6d8hPpL+nd6Lh62/aZKLB+v+OMkJSVab6uOa5DZ8XGGNV7sleXf6wf0PJvnnJM+Y6mOd4HZUr6NrJzU7tvZTr+PTjup1dO2kXsfehtOmZmuwIgAATHuTMeYVAADGhfAKAEAzhFcAAJohvAIA0AzhFQCAZsz68FpVNdX7AIyOeoV2qFcmyqwNr1X1xKr6ZJKDxun5Dqmqrqq2GY/nmwpVdVVVfXWcnuv7VXXhkJ+7qjpxPJ57Bdt9dlX9suXXgaWp16W1Wq9VtU5VfaSqvl1VD1XVrVV1XFWtORHbY/Kp16U1XK/rVtUnq+oHg3q9rqpeMhHbWhmzNrwm2SLJO5KsMU7P909J9kjy/XF6vplmjySnTuQGqmrnJJdn/F5Tpg/1OrkmpF4HPXH/kOTwJH+V5A+SnJ/kPelvIM/MoF4n10TW62fTf/vVB5K8Mv1rcHlV7TPe21sZc6Zy4zNJ13U/SfKTqd6P6arrunF5x7ksVfXEJEcneV/6b/GAEanXkU1gvf5OkpckeX3XdZ8ezPtSVS1OckJVbdV13bcnaNs0Sr2ObALrdc8k+yZ5Xdd1n0mSqvrXJM9P8ifpv3FrSkyrntdBV/hHq+rEqrq3qn5WVRdX1dOHLbdJVf31YPlHquobVXXwsGV2qaovVNVPq+oXVfWVqnrF4LF9ktw6WPRvqur7Q9bbraq+WFUPVtXCqvpcVW015PGnD7ro/7SqvlVVD1TV/17WxxpV9cyqurCqflRVDw/24aXD9rOrqvcPHltYVZ9YTtvMGxzvMVX1k6q6s6qeVVWrVdU7q+qbgy79h6vq61X1+iHr7jPYzouq6h8H7XF/VZ1VVWuP8HrsNngNPl9Va42w3Auq6trB835v6LaHHeeJw9rwwKr6dFX9vKruq6qTq+oJg9f/7sFr9/dVtf7ytj3w+iTHJJk3+JdJoF7V6yrU6+IkZyX50rD5twz+feoI6zIG6lW9rkK93pS+V/cfl8zoum5xkl8lecII6028qf6u3GHfb/v9JD9N8q0kr0nyhiQ/SHJnkrUHy2w8ZN6bk+yX5G+SdEmOGSyzTvrvL/5Skpcm+b30jb8oyTaDx18/WOdDSZ49WG+P9N/Pe1WSVw2W+ebguX57sMzTB+s9muTI9N3o2yQ5ZDB/m8Fy2yZ5YHAsByZ5RZJL0//yPnDIMXfpT4QPJvn9JLstp23mDZa7bXDMbxrM//Bgn9+ZZJ8k+yf5f4NjfeZgmX0G27lncLz7JjlusC8nD/ve4q8O/v/swWtxeZK1RnjNdhrSZn8waIe7Bvt64bDjPHFYG/4syUcH+3PqYN6tSS4aHOOxg308ZQXnzRZJ1h38/3Gvg0m9Rr1Oq3pdzn79nyS/TrLJVJ/XM3WKelWvY6jXJKsPnvtTg/VeOaXn81QX1DKKa2GSjYbMe+6g0Y8e/Hzi4IXbeti6FyR5KMlvJdl9sM5BQx5fN8lfJtlp8PM2g2UOGbLMNUn+M8maQ+atl+T+JGcMOzE+M2z7h+TxxXVh+uLaYMgyleQrg5NvtSEn3Y2jaJt5g2VfMWz++Rn8Uhky7zmDZd88rLhOHrbcvyb59vDiSrJDknvTjzNabmENafe7k8wdMu/5g+2tqLguGfL4mul/YX0vyZwh87+U5OsrcQ497nUwTdykXkdsG/U6+vPogPR/DE+b6nN6Jk/qdcS2Ua8rbqOPDJ63S3LO0NdxKqZpNWxg4J+6fnxLkqTrupvSN/jeg1kvTHJT13W3D1vvvCRPTD+m6j/Sv+BnV9XfVdX/SrJ613V/1nXdvy9ro9WPm3x+ks8nWVxVc6pqTpIHk1yZ/l3bUN9awXG8MMk/d11335Bj6ZL8XZKnpC/u0T7Xcrfbdd0buq77SFVtUFV7DD7eOXLw8PCrd68b9vMPkzxp2LynpD+h10tyeNd1j65gf/ZJ8sWu6x4ask9fSf/OfUUe25+u636Z/p3ojV3XLRqyzL1JVjRsgKmjXkemXkcwOP4FSa5N8mejXY9Vpl5Hpl6X78L058kx6XvNLxvlehNiOobXHy5j3j3p3/Fl8O+Pl7HMfw/+Xa/rul8k+d30V7W+LP27p7ur6h+qasPlbPe30neLvyP9O8+h02uy9Fisu1dwHCvcz5V4ruVut6qeU1VfSX8SXpn+D8CSC/GG32PvoWE/L87S58DT0hfwo+k/cliRDbLsgfT/vYx5wy1cxrzhF1x1o3gepo56HZl6XYbqHZ8+FF2R5GVd1z0ymnUZE/U6MvW6HF3Xfb3rumu6rjspybuSvKSqXjDa9cfbdLzbwEbLmLdJkhsH/78/yabLWGazwb/3JknXdf+V5E1VtVr68SWvTd/gC5Mcuoz1f5b+hTwtyd+u4r4PNar9HIuqenKSf07/UcwuSf6j67pfV9V2SQ4ead0RfDv9L6Q/SXJSVV3Ydd1I77B+kmUf54ZxW5PZQL2OknrtVdUa6XtbX5f+48c/HtYbxMRRr6OkXpOq2jZ9j/m5g57tJZacL5tP9D4sz3Tsef39qpq75Ieqel76MRxfGMy6Mslzq2rrYesdnP7dzFcHV/3dU1XP7rpucdd1N3Vdd2z6jzu2GCz/66Erd133YPor67bvuu7GJdNg3p+k7yZfGVcm2a+qNhg2/w3px+SMx+1gtk3/y+ivuq77Ztd1S47pZYN/V+X1vX/wEcPHk3wjyRlVte4Iy38x/Wv22DJVtWOSLVdh27RHvY6eeu2dlz64vqfrusME10mlXkdPvSa7Jjk7/UV5Qy25q8M3J2k/ljIdw+sm6W+A+4qqemOSi5PcnP4XXtK/6Pck+UJVvbmqfr+qzk5/xeGHuq5bmORrSX6ZZEFVvX5wK4vjk+ycftxG0o//SJIXVdXug/8fm2TPqrqoqv6gqvZLf2XeG5IscyzPCD6Q/mOSqwa3rHh5VX0u/RWXxwx7F7Oqbkv/jvbdVfXKqnpxVf1lkhMGjw8fbzNqgz8ob0n/ru/kERb9YPrjvKKqXl1V/zPJJel/0THzqdfRm/X1WlWvS//afz79/V1/Z9g00h9yxk69jt6sr9ckn0kfss+rqsMGbXBSkvcmOb3rultHXHsijdeVX+Mxpe8G/0z6F/OB9F3/52bIFYWD5bZI8n+T3Jfk4fTv3t44bJlt0xfm3elf6JuTvG3YMp9MP2D8pxlcOZd+QPKXBvN/nv7qwP2HrPP09B9//PGw5zokw65yT7Jj+kHNPxs837Xpx3YNXe+xqwRX0DbzBss+Ydj8vQb7+Iv0HzFclf4m4Dcn+cfBMvsM1t1v2Lp/m+THQ36+KoNbeQxroy7Ji0fYt53Sv3N/MMmP0r+T/mpWfDXk8Db8cZK/HTbvwiTfX4lzaKnXwTQxk3odsW3U67K3/en85orlZU37LW9dk3pVr5P/9zX9EIVTk9wxeK1vST92uabyfK7Bzk0L1d/M+Ktd1x041fsCjEy9QjvUKzPJdBw2AAAAyyS8AgDQjGk1bAAAAEai5xUAgGYIrwAANEN4BQCgGcIrAADNEF4BAGiG8AoAQDP+P7id7qn3RRaPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we can visualize the resulting ranks to check for uniformity:\n",
    "fig, ax = sbc_rank_plot(ranks, num_posterior_samples, plot_type=\"hist\", num_bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Result: \n",
    "The ranks appear uniformly distributed, therefore, the posterior is well-calibrated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior analysis with posterior correlations\n",
    "\n",
    "We can look for potential interactions between model parameters by\n",
    "\n",
    "- visual inspection of the posterior pairplot\n",
    "- calculating the correlation between parameter values sampled from the posterior\n",
    "- calculating conditional correlations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visual inspection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2dad9a679b81407c877616e4b7fbb59c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Drawing 10000 posterior samples', max=10000.0, style=Prog…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAJpCAYAAABGjoPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1HklEQVR4nO3de4ykV3nn8d9T1bfqnp7p8QzMeNo2QyDGeFmEI9AGJibG2uwahmCEkk2kILBQosCKOEs2a/WGbPAqCLeARFqIvAaScEmykJsQIUPAxAOBjJcEQiALGBJi2pcee+49l75Nd9XZP6q7zjnvdFVX91TVW1Xn+5FG/dZ7q1M9Vd1Pn+c855hzTgAAAP2ukHcDAAAAOoGgBwAAJIGgBwAAJIGgBwAAJIGgBwAAJIGgBwAAJGFgk+PUs6ObWN4N+InCT/ffZ6JQ9NuVcu8/T2I+X/nTXD8XffmZ6GaW+e+uN+1Ms+f1oUafCXp6AABAEjbr6QHQ7zrV60LvDnD1mu2xSahnZyvo6QEAAEkg6AEAAEkg6AEAAElgTA+AzqB6C2ivsGKLMT0boqcHAAAkgZ6eLnJo+qhm5xYlSZMTJR2buj3nFgEA0D8IerrI7NyiZqYPS5IOTh3JuTVIRph2CrU6BUVKC9ie7ESDoTCNRUprU6S3AABAEgh6AABAEkhvAanLIe1kg0PRY7dyueNtAHpGs2krKiQ3RdCTs+zgZQAA0B4EPTkLBy8DAID2YUwPAABIAj09QC/oVFl5vefJPlez59W5txXjv7fcSjONA/pMs6Xo4efNVeJbDAz6Q6t8kDZD0NOlJidK0Vw9TFYIAMDVIejpUtkAh8kKAQC4OgQ9QC/oVPnp1aawNrtu/fLLdMMDTZeiZ1JadY8xI/OmGMgMAACSQNADAACSQHoL6Cd1ZmQNZ0B25Uxqql6qarsptWZTZKF6M8lmr2eWWfSrbVZORp9nZmTeFD09AAAgCQQ9AAAgCaS3gD5iRd+97YLu7agLPNvt3WxqKVDcucNfMr+YaYP/Wyqq0tpOdztd9OhnDSYdlPnPUWF83J+2GH/e3GpQseX4vGyGnh4AAJAEgh4AAJAEgh4AAJAExvQAvabBWBu3cvmqro8XBY3PC8fqlOfO+0tGRuI2lDeePbY4sctff+FSfPBqx+5Q2o4eZAW/4KirxH0QhaFgIdHL/nNtQ0PReW51tc7NM4uZMluzJHp6AABAIgh6AABAEkhvAb2mydRNmHYKU042MhyfuLLS1HmFHWP+waV5f14x/tspnP3ZhuOu+HUDpVL0uBLcrzK/sOE1V6g0KMMP1UnZNT0zNdAmNuB/BRfGRuOD4bHR4POytBydFqa+olSXZfo0wmxXwu91enoAAEASCHoAAEASSG8B3ahRNVJwLJwZubrDH6tkK6TWBeksKZ41Obyf7YjvXT5xyjchPG9o4xRW9YbB31XB84Rd8lKcImuYggo1Wlyx3rFikwugAq0UVFIV91xT23YLwezKmc+8hemu8LOTSW9F7+kgvRVWhkmbfJYSQk8PAABIAj09OTg0fVSzc9UIf3KitMnZAACgFQh6cjA7t6iZ6cN5NwMAgKQQ9AB5qjP2JJyNVZIqwRCYeMxLPPuxC1Y8D0vOG5WBh6WytmtncCDOfhcP7Kt7LGrDQFAivujHH7gd/nkqO/fEbbjkzysEY4eyr6/pcvZ6bWvxuIawPL/hvRMuEU5SZjbkwnAw/UMwPsfGgmkgnrE7uqYSfMbcLn9eMfOzwZ2/EDxR8LOhkpmBOZiROSyVj963CczazJgeAACQBIIeAACQBNJbPWJyoqSDU0dq28embs+5RWhaowU+66gsLUWP66ZRsrOzBseKwQzK0WzKmbJttxzkzsIu9VI8I/PqLj/ofmWXb8/Kjvhvp9IJf7/BM36m5fnnTNS2h8/GJeuVEd9lP3D6rG9OqcFA/6D0vrLUYHbl4Psfpg2zM05XwtRgUEKf/b+IUgjNLPCK1uj2BTTD9mVmQ3ZBWws7x2vblQm/vXj9eHTN8oR/nw3P+ffz6MU4xVsoBTOvhyXrli1Z530rEfT0jDDIWQ9+AABA80hvAQCAJNDT04NIdfWY7VTuZFJiYdoqqrbKLvaZmUW5dv1YkCY6ez5u3g3P9PcOqqiWrtupei4d8D86Smfj11cZ8m26cLOvSLlwQ/ia4h89z/hmkKa76Yba5sp4XKkyeNGntAYf81VedimefbpcZzbqMFVVGMxUwdSrvmo0OzY6p9vSWVlh+zKZuOJeX61Y3u8/Exef7dPOyxPxZ3lhn7/J2Gx4w73ReaMnz/gmBOmtKyo7V+OZ2FNF0NODSHUBALB1pLcAAEASCHoAAEASSG8B29WoFL3F4z6uWE19TcMVl/f61ZxV8fn9cAyPJBUuBKXaq77dAxfjMQCL1/rS2OUJP8Zg/Ml47MDybv9jpbjsxzmMnHUb7peky7v8NWNzftzNYGZW2cHjZ4MHwY+vwXil97A0PVxFPnTF9y4QlrNbZlV6Sn8zwtLobh9302rZMvo14XQRUjz9w+q4f68u7wymRsj8OBkKht6NzPnPWOl4PF4tLId3l/37sTA6qvhEf49w7E9q6OkBAABJIOgBAABJIL0FbEWdBUIbntdAuHhoNAtwJiUTlmAPXH/AHxiMZwt2Fy/6B9f6NJYr+XRP8clT4SWq7JmobS8926fRBubjLvBL1/q27prxr/34oTi1dHm3P1Z62v+IWR3z3fCX98X3Hn00mHG64NswOpt5fWHaKUhP2XDcBq3486LS9HAW58z3OPq/CBd1HBmJznNzPu8QzZTdKNXV7PumF0Wl2l0+a/LVqpPOkqRCMF2EjcWppfIu/3jhmUF6a7e/n/3Yueiai6d9imzknP8c2Urm/RPOMj7k722Z9JZbDtK5Cack6ekBAABJIOgBAABJIL0FbEWzqYk651mmyii6JLuwZShIj5SPn/C7x+IFOQvBzK/uhJ+ptfzcICV23TPia+Y3rmJaHYt/PKwGT/X0v/N/L5XH4tc68rRv68rOoGJr0XepjzwRfx9c8OfX+Hd9+mj+ufGs0APBDM3z+4Nu/XLcRT/2tE9dDT8epA3mgvTfqTjNF97BzftFHbNpsGYqw67Qbymtevo9VZJ9feHis4u+ClKZ6q2VXT5Feuk6/2aff4H/7A1djj9vI4/7z0hx2b9/Tty6Jzpv/+eDtOolv8BvlOrWNiu2ui0N1oL2EPT0uHBJivXHLEsBAMCVCHp6XDbAYVkKAAA2RtADtFvYBZ5d1DJMezRZ8RWyUpzeCrvYXVC9NXDWd3tbZjHO8gHfXT580qd1Fg/EXfQWzEE4NOe7mQun4nYXgl704XP+vNGffLq2/dQ/xym26553srZ95sn9/l4rcRf28UP+uUongnufjE5T4XI8YeK6bJd/KKzeCtNWhUwlTjSpYavTVv1c5dXnbNhXDLrdcVq2uOQ/FJd3BpWFZ3yqdP+BM9E1jz/HHzsVVC3u+Vb8mXCn/YSdFqR/GqZem00TdUNKK9SC9jCQGQAAJIGgBwAAJIGgBwAAJIExPcB2Zcfg1BuD0WBsRjSjb3a8TxOiMllJFs4Ke9GP41HB/33jdsRjVFYmfDltccGPPSiPxH8TjZ7w42ROv8qPa6msxrPUlr7r7zcYvKRzf+vH6gz+23hc0Ym5cd++m/z+nd+P773/73wbzv+Q//6PPx6X3ZdL/tjQJT9OyXZP+O3M9y6c9Tr6P8ssOBpOE1C+EJzXijE4vTiOp9vKmtspMyOzBYvehrMhZ89b3Oc/E+Uf9u/H8oK/frmc+Xky5+8Xjl/b8URmaotwLFEw63I4Rk2SXPje6sf/pwazZYfo6QEAAEkg6AEAAEkgvQVsV7OpiAZlyHUXqWxQvh6mwVw5Ls125y/Utu3APr89H6RySsPhJVFK6/Ju36WeLRdfmvRtGvi+T/Gs7IrbUA7W51ze64/d+COP17Y/87zPRNf83nmf+nrnudfUtgfn47/LVkaDx8HTXjgYv6aovP6cL4+3RZ+qsuU4zVcMUo2VYFHRbOnvQDDrdfi/VM5MBdCTqart6MdUSbPClHQhKBcPFviVpPl9wQzmK377JTf9oLb9Jz/0UHTNjd98S217xWettbIj/rU9NO6nlrAFnzoLU9rSJgvi9oMm34f09AAAgCQQ9AAAgCSQ3gLabTtpjuw1ddJd0ezAkorXXetvUfLpmoXn7K5tD52P0zXF+aDbO0hvucxTjs/6Nq3s8Aft5vnovMp8UEEWVHZ9/++eVdt+/fBt0TVv3e+79otBSmt5Iq7IWBn3jy+P++7socxEy3PP98eWrvGVYeNP+tew47GF6JpCsBhpmELMVsFUTgcz5w76NEZ25uZKsGhpMqmuftcghVJ5tl/Ud3F//F4oj/j37YF9c7Xtb5/wad1vTGYW/g1StGGa+MKz4l/bY/8YVyHWmrq08ULCqSPo6TPhAqQsPgoAgEfQ02fCIIfFRwEA8BjTAwAAkkBPD9CNMmN4suNK6qkEKy6Hf9GMFINy2qH4Y18e8+N4xv7peG176Xn7o/OWd/vrKsHkswN/Nx6dV3mpL/deOOPHNgyd8tf/w+z10TX/ff51te3R5/jrl87tVj0DC/41WTkeazFyshAc8/vn9/nvY3E5XqF+dMWXtheW/Ky32bERleCYGq1kjb4XzcK8Un/c1uroxvuvm/Dv9enjr4wPBh/gwQv+QaFcf1xROCPzdmZ4TwE9PQAAIAkEPQAAIAmktzrg0PRRzc75ssLJiVKDs9HVGsyu3ErNprOuuG4omAm24stciyfmatvhgpnVg8FipMEinCOPno5Ocz/s0z+lIH106VnRaao8urO2PRykoIZ8T74WZnaEl+gJ+cfl3T5ltP9f49mez97s77f/7/33/+zz4x9lo0/7FEBY9h7O1Jydcdrmg5TWSoO0VZ33QGFkJDotO5Mzeki9RVSzaedR/1mqDPjPROlEZjHbEf/+Pv50kLIt+Hu7pfjeY3O+DeVgwvFrvpO595lzQYOCfgy3KlyJoKcDZucWNTN9OO9mAACQNNJbAAAgCfT0AHmqkyrJVl5Es/0GqZeokkiSC2YBdkF6JaxAKgYLkUpSZbevvnIH/YzOq+PxIp5DZ3y3+u6ggmTnY/HfTsXlOCW17vIu/+NmZTzuyg9TX8VHw8Ua4xTUrn/d8Na65pG4K3/0B37h1cUbgtcXVLEt74rbMLrov0eV8savQYpTj2HrrkhnMQtzfwhSXYWR+DPhlv1s5oUFv71wcCI6b/cj/p1y6YK/RznINA+fU2Ti+/79tLg3WOz35IXovEqw0Km73OeLirYAPT0AACAJBD0AACAJBD0AACAJjOkBtiIcp5Fd+byVYzgy9wrH5DSaadUFY1HcZT++JxwT5C7Fq6K7vb7EPJzZePgHp+LzFvz4oeELfpzM0rOvic4bnIvHGdWuedK/hsFLE9Gx+f1+HM9qMM6hdCYeq+OC8QuLe/z3f+ejcRnv5X1jte1wlunytb6to4uZMTjB97Ww05cYVy5citsQfv8bvB/CEvZovA9jfXpWdsxMOCNzeZd/447OzMXn3bCrtj18wb+HLx3wv4ILmeE4xWX/Ppn4bvD+nrtYt0024D9HbjXz/m6wQnxK6OkBAABJIOgBAABJIL0FbFc2TbGd2ZpbnOqwoKQ2LF+vBNuKs1sqlsNS+SA9Fs7unLmfm/M15kMzT8b3u86XvYcLoNoun0YbOhOX/l7e6b93e76Wqd0N2Lzv5h884FNVA2fjFzU4E6SkghLz4sk5/xqyab5rJvzzLPprCpm0lQumCYi+r5n/ywrVw/0hnOU4M1N6mOos/ov/HLjM9AWjF/x7bekm//nY+48+VRWmbiVp4Cn/OXAjPo1mmc9l3ZQW6awN0dMDAACSQNADAACSQHoLaJVtpKps0Hdb160KkuRWmsuV1Eu31K0kklSZjyufas8ZVIxJcerMgh72bHe7u+C77C1czDSYvbZ4Nq6I2vn4Sf/gGcGCjCcyi54G36OB8/HMtNF5QRoirLAJ26AgnSVJFs7IfPqMPzAYv75wRuxG/+fRzM0NKu7QO7LVW4XhIJ286NOe4UKkWcOPBSnfS0HKODsLePh5W/LPWwlSy9XrWvjesjjF1o8pMnp6AABAEgh6AABAEkhvATmqm7ZqNPFh9lg9hTC94rvOC5l0VJjusmL9e1cu+rRVccJPthZWdUnxIqjheeUzvlu/2CBd5x593Lcnu8Bj+Dp2+AkIK+fm6t8veH3hNcqk+VaPP+3bt9u3Ozs5YaRBxV70fxuet50qP3RWkNaxYpjyyXw+CkG/QcW/NysX4/eMBWlVW/DpZDfgfwW7hfhzZAPB5zf8jGXTWfV+NrhtvLf6MJ2VRU8PAABIAkEPAABIAkEPAABIAmN6gG7UihmdwzE9wfgSt9LceY3GDpUzZbOhsAw/Oi+4X7PXK1PGG5bkh+Xw2bLy8FhYku/CEvrMOKBwrFM0+3RmDEU07qnJBWij8vUmpx9Ad2i4wG84rUOD90JY6l4oBeXsi/69GU6tIEmVC/HCorX9y5mpJMJxQaur2dORQU8PAABIAj09iTg0fVSzc9W/KiYnSjo2dXvOLQIAoLMIetokG2TkYXKipINTR2rbM9OHJam2D9vQDeXGTaa0tpUia8Hrq5u+aTYVFy6AmimHD5WDUvIrSu3D2WzrpboapC3CVFejEv9mpw8gpdWFwtmHmy3VvuK8yoanNVqMuFIvpZVJTbmV4HGwGGn2vU5Ka2sIetpkdm6xFmTkhd4cAAA8xvQAAIAk0NMDbEW3z6DbTPsazfac0+urtzjnFamlYLbnkKtkK6w2ruxyS0FFTIPUVDgTdDkzI3O9WatJYfWY7aS0trkgpxXqzOocpnIrmXsVMs9V7zxsCT09AAAgCQQ9AAAgCQQ9AAAgCYzpAfpVvfLzLhyXFJWPB+2rLNUv/W00y3FUVl5nHFCjsuJGM0Znxw9tdP2G90dv2k5pe0b4/q43fk0uU3pe73mz44qwJfT0AACAJBD0AACAJJDeAnpZF5afb8s2Zo+um2ZSJoXQisVbt3N9N8zeje1pQUorEtwjmkE5fI+4zHukFc+LK9DTAwAAkkDQAwAAkkB6C+hlrV58tJc0eH3MjoyuFabO+vFz2eXo6QEAAEkg6AEAAEkg6AEAAElgTA/Qr/p9vEC3v75ubx/qa2e5+NXem1L2q0JPDwAASAJBDwAASALpLQD9rd9L99G7Wj3zMzZFTw8AAEgCQQ8AAEgC6S0A/Y2UFroVKa2OI+hJ0ORESQenjkSPj03dnmOLAABoP4KeBGUDnDAAAgCgXzGmBwAAJIGgBwAAJIGgBwAAJIExPS10aPqoZucWJVUHBwMAgO5B0NNCs3OLmpk+nHczAADABkhvAQCAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJDAjMwC0QqEYP66U82kH+oOZ33Yuv3b0GXp6AABAEujpgSYnSjo4daS2fWzq9pxbBABA6xH0XIVwVXWpd1dWD4Oc9eAH6HutTkeRzkIzmk1bkdJqC4Keq8Cq6gAA9A7G9AAAgCTQ0wMgTaSjkAfSVrmipwcAACSBoAcAACSBoAcAACSBoAcAACSBgcxbFM7N06vz8jTCRIUAgH5F0LNF/T43DxMVAgD6FektAACQBHp6UBepLgBAPyHoQV2kugAA/YSgBwCywsVIMzM32+BQbdutXL7q+yFdNuB/BbvV1czBJhcmxZYwpgcAACSBnp4m9HuZOgAAKSDoaUK/l6kDyGiQgmo6pdXk/ZCuK1Ja0UFSWu1AegsAACSBnp46SGkBANBfCHrWhEGOVA10SGl54Zw964+ZtwcA0EsIetYwbqexbIBzaPooExcCeQlL4JGWsJRdYuzPFhH0YFuanbgwmyYkOAIA5GXbQU+v/DLLpq3qYdzO9mVTX9lj6z1ozOoMAMiTObrGAABAAihZBwAASSDoAQAASSDoAQAASSDoAQAASWhYvWVmn5W0t0NtabW9kk7n3YgOSOl1ftc5d0feDQEA9Ka+rd4ys685516cdzvajdcJAEBzSG8BAIAkEPQAAIAk9HPQ88G8G9AhvE4AAJrQt2N6AAAAQv3c0wMAAFDTV0GPmRXN7D4ze8rMLpnZn5nZvrzb1Wpmts/MPrr2OufM7HNm9oK829VOZvajZrZqZrfl3RYAQG/qq6BH0r2S3ijpDZJeLuk6SX+eZ4NazcwKkj4p6UZJd0p6maTzkh4ysz15tq1dzGxM0h9IKubdFgBA7+qbMT1mNqTqJH13O+c+srbvoKQfSDrknHs4v9a1jpndIunrkm52zj2ytm9Y0llJb3HOfSzP9rWDmX1A1SDvNkmvcM59MdcGAQB6Uj/19LxI0rikL67vcM7NSJqRdGseDWqTxyW9WtL3gn0VSSZpdy4taiMze5Wkw5LuzrstAIDe1nAZih5z3drX2cz+45Ku73Bb2sY5d0bSkczuuyWNSHqw8y1qHzPbK+l3Jb1J0rmcmwMA6HH91NMzKqninFvJ7F9WNSDoS2b2Gkn3Sfrt9XRXH/mApE875z6bd0MAAL2vn3p6FiUVzGzAObca7B+WNJ9Tm9rKzO6S9CFJn5B0T76taS0ze6OkWyS9MO+2AAD6Qz/19Dyx9vXazP4DujLl1fPM7O2SPizpAUlvcM5Vcm5Sq92lasryaTO7JD+G6a/M7IHcWgUA6Fn91NPzTUkXJf24pD+UatVbByV9KbdWtYGZ3SPpnZJ+wzn3m3m3p01eL6kUPN4v6cuSfl7S53NpEQCgp/VN0OOcWzaz+yW918xOSzop6X5Jf+Oc+0q+rWsdM3uhpHdJ+n1JHzKz/cHhi865vkjlOeei3jkzW1rbnHXOncyhSQCAHtdP6S1J+nVJf6RqT88XJD0m6adybVHr/ayqk/S9SdJTmX9vy7FdAAB0tb6ZnBAAAKCRfuvpAQAA2BBBDwAASAJBDwAASAJBDwAASAJBDwAASAJBDwAASAJBT5PM7DYzc2Z2x9rju9Ye39SB536PmfXNBIsAAOSBoGf7jkh6qaSZdj7J2pITv9rO5wAAIAV9swxFpznnTkk61a77m9mNkn5L0islnW/X8wAAkAqCng2YmUn6b5LerOqq7Q9L+mjmnLtUXeX8+c6575rZvaquDP6fVV0b6yZJP5D0XyV9X9LvSPoxVQOl9zjnfmeTZnxQ0i5Jt6m6uOjIVb8wtAJTmKObWN4N+InCT/OZQFf5fOVP634uSG9t7F2S7pP0CUl3SvqGpA81cd0zJT2gag/NayWVJX1c0mdUXRn8dZK+I+n9ZvaSTe51t3PuFufc326j/QAAIIOengwz2ynpVyTd75z7tbXdD5rZuKRf2OTykqRfds59cu1e+yR9RNL7nHPvWdv3bUmPS3qZpK/Wu5Fz7p+u5nUAAIAYPT1XeqmkIUmfyuz/eJPXHwu2T6x9DSuvTq993b31pgEAgO0i6LnSnrWv2UHKTzV5/YUN9s0H2+S/AQDIAUHPldaDnf2Z/Xs73RAAANA6BD1XeljSgqSfyey/M4e2AACAFmEgc4Zzbt7M3iHp3WZ2QdVJCG+V9Iv5tgwAAFwNgp4NOOfea2YXVa3ierOkb6paufWJXBsG9IlD00c1O7coSZqcKOnY1O05twhACsw5xtUCW8AHpgUOTh3RzPThK7axZUxOCGQ0mpyQnh4AuZqcKOng1JHoMT0/ANqBoAdArrIBThgAAUArUb0FAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSwDw9AAD0KgsmH7agH6NS7nxbegA9PQAAIAkEPQAAIAmktwAA6BWWWUszSGlZwR9zLnMei4tLoqcHAAAkgqAHAAAkgfQWAAC9KqjSsqGRuqe51dVOtKbr0dMDAACSQNADAACSQNADAACSwJgeAAB6hBWL0WNXLgfbFX/e0FB8XjimJyx7T6yUnZ4eAACQBIIeAACQBNJbAADkodk0U8GntLJpK5XLGx6zzMzNNuB/3adcvk5PDwAASAJBDwAASALpLQAA8tBk5VRYsWWlzKzLi0t+O6zkyj5VJa0qrXoIegB0xKHpo5qdW5QkTU6Ucm4NgBQR9ADoiNm5Rc1MH867GQASxpgeAACQBHp6AADoMmGJeWFiV93zwhmZwzE9hR1j8YlLyy1rWy+jpwcAACSBoAcAACSB9BYAAF0mLDG3kWG/f0dc+Vgc9SXsq08cr21X5hfjG1bKAj09AAAgEfT0ALgq2fl3jk3dnnOLAGBjBD0Arko4/87BqSM5twboD8VdO2vbK9fvqW3PT8YzMk/83ydr2zbof6VnFxx1q00ubtrnSG8BAIAkEPQAAIAkEPQAAIAkMKYHQMtMTpSicT0MbAaalBmDY2Ojte1zz/PbA0uZ8TjBLMzhOB63uhqfl/A4nhBBD4CWyQY4DGwG0E1IbwEAgCTQ0wMAQM4KpXimZRfMtLy826etRv8lM7NysVjbrCwHi4qSztoQPT0AACAJBD0AACAJpLcAbFl26Yl6wmquRucBKSqM+qosK8UzLS8855ra9sq43794TTE6b2x+oT2N61MEPQC2LFx6ohHK1QF0E4IeAF0l2ztE4ASgVQh6AHSVMMhhnh/0s8riYm17YP8zo2MXDvpfz8Nn/X7LFmUFkxNSsbU5BjIDAIAkEPQAAIAkEPQAAIAkMKYHAIBOKfiS8+LOHbXtCy/aF512eZefhXnh2kpt+9l/sRKdd8XComiInh4AAJAEgh4AAJAE0lsAAHSKq2y4e+658UzLLuiSGHvCPxicW4ovrGx8P2yMnh4AAJAEgh4AAJAE0lsAAHRIuMioOzhZ217NrMc7dN5vjz3tU1i2GqezXJn01lbQ0wMAAJJA0AMAAJJA0AMAAJLAmB4AANrJ/OzKNuh/7R5/+URte8dLTkeXnP2Xa2rb13zXr6RuT5+JzquEq6xjU/T0AACAJBD0AACAJJDeAgCgjazoZ1su3/Ss2vaF5/rU1KM/8ifRNf9m+edq28sT47Xt0vhYfPMTJ1vVzCTQ0wMAAJJA0AMAAJJAegsAgDayAf+r9vwP+RmZh/bP17ZPluejaxaP76ht750PZmS+tFD33q7i/IEKVV0boacHAAAkgaAHAAAkgfQWgJ50aPqoZucWJUmTEyUdm7o95xYB6HYEPQB6QhjkSNVAZ2b6sCTp4NSRvJoFbMqGhmrbA0t+3M3goB93c+e33hBdM3LSl7mPnLpU267Mx2N6XDgjs3NCYwQ9ALrW5ESpFtCEQQ4AbAdBD4CmZNNJnUDKCkArEfQAaMrs3CI9LcA22E4/o/LpF/q01eVLw7Xt5e/tiq458PWV2vbAifO17dVMeouU1tZQvQUAAJJA0AMAAJJAegsAgBaywaHo8cr1e2vb7uaL/sBxv3ioZSZQHpj3Oyqnz7a2gQmjpwcAACSBoAcAACSB9BYAAC1kI8PR45Mv9mmsYtFXYq3sXa5t7/5qfE3hsk9vucuX/QEWEr0q9PQAAIAkEPQAAIAkEPQAAIAkMKYHAICrZAP+12lhx1h0rBwM1ykWK7Xt4Uf8ci6uEM+sPPDPT9a2K0Kr0NMDAACSQNADAACSQHoLQM+bnCjp4NSR2jars6PTXMWnp9x4nN4ae9onqM58d6K2vfOUv2bPP5yJb7i66u8XlqzjqhD0AOh5YZCzHvwAQBbpLQAAkAR6egD0FVJdyENhbLS2Xb4mTm8Vg+zUnv/nU1rjj/sZmXV5Jbqmcmm+tQ2EJIIeAHUcmj6q2bnF2uPJiVKDs7sHqS4A9RD0ANjQ7NyiZqYP590MAGgZxvQAAIAk0NMDAMB2FIp+c+d4bdsuLkenFS+P1LaH53z5+uCZYNzOybhk3QUl62gdenoAAEASCHoAAEASSG8BALANNuh/hVb27PTbo0PReQMLPqU1dMZXRNr5S7Xt8sJCO5qIDHp6AABAEgh6AABAEkhvAQBwlWzRT7tc2R1P5Dkyc84/KPq+Bjfvq7fCBUvRPvT0AACAJBD0AACAJBD0AACAJDCmB0DfCldcX3/MqutoFSv6GZk14LcHzsczMpcn/ArsxTlfml5ZXGpf47Ahgh4ANeHK6r2yqnoj2QCHVdeBtBH0AKhhZXUA/YygB0AywnQXqS5smVn9YwU/RNaWV6JDxZWyP3YxWGS07PerEmyjbQh6gMT1W0qrkTDIIdUFpIegB0gcKS0AqSDoAQCgCTY0lNnh0112IUhbra5Gp7nLwWzNl33qy5VJaXUa8/QAAIAkEPQAAIAkkN4CAKAJbjmedFDBIqHu4kW/fyD+1eqWguuc23gbHUFPDwAASAJBDwAASAJBDwAASAJjegAkidmZ0ZSwLD1cYFSSFX2/QeWSL1kv7NoZ3yMYu1NZWBDyQ9ADJCacgVnq/1mY62F2ZiA9BD1AArJLTTADM4AUEfQACWCpCaAFMumtkAvK18NUlyS5ldXs6cgJA5kBAEASCHoAAEASSG8BAFBPUHmVTVNF8ym7it8MFhiV4qqv4DTkgKAH6FPZwcsAkDqCHqBPMXgZAGKM6QEAAEmgpwcAgGZUynUPWbCyulvNjP1ZpWS9W9DTAwAAkkDQAwAAkkB6CwCAegrBLMzZ9FYhLEUPCtiDRUqrB53QHejpAQAASSDoAQAASSC9BQBAPQ0qtupOr0w6q2sR9ABI3uRESQenjkSPj03dnmOLALQDQQ/QR1h6YnuyAU4YAAHoHwQ9QB9h6QkAqI+gBwCA7XANytTRlajeAgAASSDoAQAASSC9BQAZYTUXlVxoCmXqPYGgBwAywiCHSi6gfxD0AD0sLFGXKFMHgEYIeoAeRol6+5HqAvoHQQ/QY5iAsLPCIOfQ9FECIKCHEfQAPYbenfww1gfobZSsAwCAJNDTA/QAUloAcPUIeoAeQEoLAK4eQQ8AbENY1bX+mIHNQHcj6AG6EPPvdL9sgENlVx8JFw9lpuW+QtAD5Cg7Vmf9FyXprN5Tr7Q9i4AIyA9BD7AF2R6YesJfbI2umZwo1YKbbE8BelejoIZSdyA/5ui6AwAACWCeHgAAkASCHgAAkASCHgAAkASCHgAAkASqt4AtMLPPStqb09PvlXQ6p+fuVql/T0475+7IuxFAr6B6C+gRZvY159yL825HN+F7AmArSG8BAIAkEPQAAIAkEPQAveODeTegC/E9AdA0xvQAAIAk0NMDAACSQNADdDEzK5rZfWb2lJldMrM/M7N9ebcrT2a2z8w+uvY9mTOzz5nZC/JuF4DuR9ADdLd7Jb1R0hskvVzSdZL+PM8G5cnMCpI+KelGSXdKepmk85IeMrM9ebYNQPdjTA/QpcxsSNWJ9+52zn1kbd9BST+QdMg593B+rcuHmd0i6euSbnbOPbK2b1jSWUlvcc59LM/2Aehu9PQA3etFksYlfXF9h3NuRtKMpFvzaFAXeFzSqyV9L9hXkWSSdufSIgA9g2UogO513drX2cz+45Ku73BbuoJz7oykI5ndd0sakfRg51sEoJfQ0wN0r1FJFefcSmb/sqq/5JNnZq+RdJ+k315PdwFAPQQ9QPdalFQws2yP7LCk+Rza01XM7C5VB3X/saR78m0NgF5A0AN0ryfWvl6b2X9AV6a8kmJmb5f0YUkPSHqDc66Sc5MA9ACCHqB7fVPSRUk/vr5jrXrroKQv5dOk/JnZPZLeKek3nHO/5ChBBdAkStaBLmZm05LuWvt3UtL9kpacc7fl16r8mNkLVS1Z/6ikt2cOX3TOJZ/2A1AfPT1Ad/t1SX8k6Q8lfUHSY5J+KtcW5etnJRUlvUnSU5l/b8uxXQB6AD09AAAgCfT0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0ANgWM7vNzJyZ3bH2+K61xze16fluMbO/NLMTZnbOzL5kZre347kA9CeCHgCtckTSSyXNtPrGZvY8SV+WNCHpzZJ+TtKcpL82s59s9fMB6E/Z1ZsBYFucc6cknWrT7X9F0oKkO5xzlyTJzD4n6Z8kTUn6dJueF0AfoacHwKas6h4ze9TMFs3sIUk3ZM6J0ltmdq+ZzZjZq8zsG2a2ZGaPrD2+0cweNLMFM3vMzN66SRO+Jem96wGPJDnnypL+WdJki18utsbxj39d9K8hgh4AzXiXpPskfULSnZK+IelDTVz3TEkPSPotSa+VVJb0cUmfkfR5Sa+T9B1J7zezl9S7iXPu/c65d4f7zOwaVVeg/9bWXgqAVJHeAtCQme1UNb10v3Pu19Z2P2hm45J+YZPLS5J+2Tn3ybV77ZP0EUnvc869Z23ftyU9Lullkr7aZJsGJP2+pJ2qBmQAsCl6egBs5qWShiR9KrP/401efyzYPrH29SvBvtNrX3c3czMzK0n6U1V7nO52zj3cZDsAJI6eHgCb2bP2NTtI+akmr7+wwb75YHvTPPw6MzugavB1i6Q3O+c+0Oy1AEBPD4DNrAc7+zP793ayEWb2Akl/L+n5kl5LwANgqwh6AGzmYVXLxX8ms//OTjXAzG6Q9JCkQUkvd879ZaeeG0D/IL0FoCHn3LyZvUPSu83sgqqTEN4q6Rc72Iz/rWol2C9JGjKzHw2OXXbOfb2DbQFa7tD0Uc3OLUqSJidKOjbFZOPtQNADYFPOufea2UVVq7jeLOmbqlZufaLdz21mY5Jeufbw/RucckJXpt6AnjI7t6iZ6cOSpINTR3JuTf8i6AHQlLUxNNlxNH8cHP+IquXo64/vlXRv5h6flWSZfUvZfZnj8yIVD6AF+EECAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSMJB3AwAASNGh6aOanVuUJE1OlHJuTRoIegAAyMHs3KJmpg/n3YykkN4CAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJIOgBAABJYBkKAAC6yORESQenjkSPj03dnmOL+gdBDwAAXSQb4IQBEK4O6S0AAJAEgh4AAJAEgh4AAJAExvQAANAhh6aPanZuUVJ1gDI6i6AHAIAOmZ1b1Mz04bybkSzSWwAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAkEPQAAIAksQwEAQBebnCjp4NSR2vaxqdtzblHvIugBAKCLhUHOevCD7SG9BQAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkkDQAwAAkjCQdwMAAOhXh6aPanZusfZ4cqKUY2tA0AMAQJvMzi1qZvpw3s3AGtJbAAAgCQQ9AAAgCaS3AADoEZMTJR2cOlLbPjZ1e84t6i0EPQAA9IgwyFkPftA80lsAACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJBD0AACAJA3k3AACAfnJo+qhm5xYlSZMTpZxbgxBBDwAALTQ7t6iZ6cN5NwMbIL0FAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSQNADAACSwIKjAABcJVZW7w0EPQAAXCVWVu8NpLcAAEASCHoAAEASCHoAAEASGNMDAEAPmpwo6eDUkejxsanbc2xR9yPoAQCgB2UDnDAAwsZIbwEAgCQQ9AAAgCQQ9AAAgCQQ9AAAgCQQ9AAAgCQQ9AAAgCQQ9AAAgCQwTw8AAFsUrqoudcfK6tnJCsP9TFpYRdADAMAWdeOq6vUCGyYt9EhvAQCAJBD0AACAJJDeAgCgCeE4nm4Yw4OtI+gBACCQHaS8bnKi1HXjeJpRb4Dz+rGUBjmbcy7vNgAAALQdY3oAAEASCHoAAEASCHoAAEASCHoAAEASqN4CAGybmX1W0t6829HAXkmn825EB6X2eqX4NZ92zt1R70SqtwAAfcvMvuace3He7eiU1F6vtLXXTHoLAAAkgaAHAAAkgaAHANDPPph3AzostdcrbeE1M6YHAAAkgZ4eAACQBIIeAEBfMbOimd1nZk+Z2SUz+zMz25d3u9rFzPaZ2UfXXu+cmX3OzF6Qd7s6xcx+1MxWzey2zc4l6AEA9Jt7Jb1R0hskvVzSdZL+PM8GtYuZFSR9UtKNku6U9DJJ5yU9ZGZ78mxbJ5jZmKQ/kFRs6nzG9AAA+oWZDak6Ud3dzrmPrO07KOkHkg455x7Or3WtZ2a3SPq6pJudc4+s7RuWdFbSW5xzH8uzfe1mZh9QNeC7TdIrnHNfbHQ+PT0AgH7yIknjkr64vsM5NyNpRtKteTSozR6X9GpJ3wv2VSSZpN25tKhDzOxVkg5LurvZa1iGAgDQT65b+zqb2X9c0vUdbkvbOefOSDqS2X23pBFJD3a+RZ1hZnsl/a6kN0k61+x19PQAAPrJqKSKc24ls39Z1UCgr5nZayTdJ+m319NdfeoDkj7tnPvsVi6ipwcA0E8WJRXMbMA5txrsH5Y0n1ObOsLM7pL0IUmfkHRPvq1pHzN7o6RbJL1wq9fS0wMA6CdPrH29NrP/gK5MefUNM3u7pA9LekDSG5xzlZyb1E53qZrGfNrMLsmPZ/orM3ug0YX09AAA+sk3JV2U9OOS/lCqVW8dlPSl3FrVRmZ2j6R3SvoN59xv5t2eDni9pFLweL+kL0v6eUmfb3QhQQ8AoG8455bN7H5J7zWz05JOSrpf0t84576Sb+taz8xeKOldkn5f0ofMbH9w+KJzru9Ses65qMfOzJbWNmedcycbXUt6CwDQb35d0h+p2tPzBUmPSfqpXFvUPj+r6sR8b5L0VObf23JsV1dickIAAJAEenoAAEASCHoAAEASCHoAAEASCHoAAEASCHoAAEASCHoAAEASCHoAAGgBM7vNzJyZ3bH2+K61xze16fleYmZ/bWbnzeyUmX3YzPa147n6BUEPAADtcUTSSyXNtPrGZnazpKOSnKoTFL5V0q2SvmBmQ61+vn7BMhQAALSBc+6UpFNtuv2vqrqi/J3OuQVJMrOTqgZCr5T0qTY9b0+jpwcAgC2yqnvM7FEzWzSzhyTdkDknSm+Z2b1mNmNmrzKzb5jZkpk9svb4RjN70MwWzOwxM3vrJk34H5L+w3rAs2Z57etIy15onyHoAQBg694l6T5Jn5B0p6RvSPpQE9c9U9IDkn5L0msllSV9XNJnVF0h/HWSviPp/Wb2kno3cc7NOue+IUlmVjKzQ6ourDqjaloNGyC9BQDAFpjZTkm/Iul+59yvre1+0MzGJf3CJpeXJP2yc+6Ta/faJ+kjkt7nnHvP2r5vS3pc0sskfbWJJs2oGkwtSfpp59ylLb2ghNDTAwDA1rxU0pCuHDfz8SavPxZsn1j7+pVg3+m1r7ubvN8bJd0h6S8lfcrMXt/kdckh6AEAYGv2rH3NDlJ+qsnrL2ywbz7YdltpjHPus865z0n6T5L+XtI7tnJ9Sgh6AADYmvVgZ39m/95ONcDMXm1mrwj3OeecpK9Lur5T7eg1BD0AAGzNw5IWJP1MZv+dHWzD2yT9npkNr+9Ym5/nFaoOqsYGGMgMAMAWOOfmzewdkt5tZhdUrZa6VdIvdrAZ/1PSQ5L+wsz+l6RhSf9F0nMk/fsOtqOn0NMDAMAWOefeK+ktqk4E+GlJ/1GbV2618vm/pGqvzqCk/yPpY6pOVvhjzrkvd6odvcaqKUAAAID+Rk8PAABIAkEPAABIAkEPAABIAkEPAABIAkEPAABIAkEPAABIAkEPAABIAkEPAABIAkEPAABIwv8HE2nt4t6jI0gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "posterior_samples = posterior.sample((10000,))\n",
    "pairplot(posterior_samples);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can recognize a potential correlation by looking for elongated shapes in the 2-D posterior marginals. \n",
    "\n",
    "The posterior 1-D marginals are relatively broad, similar to the prior. However, in the 2-D marginals we see strong signs of correlations, indicating parameter interactions. \n",
    "\n",
    "We can quantify these correlations by calucalting the Pearson correlation coefficients between posterior samples: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -0.29635301,  0.25109896],\n",
       "       [-0.29635301,  1.        ,  0.81718692],\n",
       "       [ 0.25109896,  0.81718692,  1.        ]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.corrcoef(posterior_samples, rowvar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The correlation matrix confirms the strong correlations between all three parameters, a negative correlation between `dim 1` and `dim 2` and a positive correlation between `dim 2` and `dim 3`. \n",
    "\n",
    "In other scenarios, the correlation might not be visible as obviously as here. \n",
    "For example, in a high-dimensional posterior, correlations might not be visible in the 2-D marginals shown in the pairplot. In that case we might be able to still seem them when calculating conditional correlations, e.g., by calculating the correlations between posterior samples that were sampled by conditioning on other parameter values. See the corresponding tuturial in the `sbi` package for details: \n",
    "https://github.com/mackelab/sbi/blob/main/tutorials/07_conditional_distributions.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The end\n",
    "\n",
    "This ends the workflow for applying SBI to wiring rules in the barrel cortex. \n",
    "\n",
    "The workflow of `simulation`, `inference`, and `analysis` can be applied to any other computational connectomics model as well. All you need is\n",
    "\n",
    "- a model with free parameters that can simulate data given parameter\n",
    "- corresponding observed data\n",
    "- a prior distribution over the free parameters, e.g., a uniform prior in the ranges of parameters values you find plausible\n",
    "\n",
    "Just install the `sbi` package from https://github.com/mackelab/sbi and follow these three tutorials to set up your own workflow. \n",
    "\n",
    "Do not hesitate to create an issue in this or in the https://github.com/mackelab/sbi repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}