{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting a model's parameters with run-at-a-time optimization\n", "\n", "In this notebook, we'll fit a simple compartmental model to disease propagation data. We'll use a standard optimizer built into the python `scipy` library to set two independent parameters to minimize the sum of squared errors between the model's timeseries output and data from the World Health Organization.\n", "\n", "## About this technique\n", "A run-at-a-time optimization runs the simulation forward from a single starting point, and so only requires an *a-priori* full state estimate for this initial condition. This makes it especially appropriate when we only have partial information about the state of the system. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ingredients:\n", "\n", "We'll use the familiar `pandas` library along with `pysd`, and introduce the optimization functionality provided by `scipy.optimize`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import pandas as pd\n", "import pysd\n", "import scipy.optimize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model that we'll try to fit is simple 'Susceptible-Infectious' model. This model assumes that everyone is either susceptible, or infectious. It assumes that there is no recovery, or death; and doesn't account for changes in behavior due to the presence of the disease. But it is super simple, and Çso we'll accept those limitations for now, until we've seen it's fit to the data.\n", "\n", "\"Stock" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. image:: ../../../source/models/Epidemic/SI_Model.png\n", " :width: 600 px" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll hold **infectivity** constant, and try to infer values for the **total population** and the **contact frequency**." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "model = pysd.read_vensim('../../models/Epidemic/SI_Model.mdl')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll fit our model to data from the WHO patient database for Sierra Leone. We see the standard *S-Shaped* growth in the cumulative infections curve. As our model has no structure for representing recovery or death, we will compare this directly to the **Population Infected with Ebola**. We format this dataset in the notebook [Ebola Data Loader](https://pysd-cookbook.readthedocs.io/en/latest/data/Ebola/Ebola_Data_Loader.html)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgwUlEQVR4nO3dd3wUdf7H8dem7KZXQkIJEDooIJ2IUhQJiAjYlTtRUU8FESwov1Owg6CnIip6d2JFsOGpCIgIKIiAIIiCNGkKoaf33e/vj0kWIjUQmGzyfj4e+9jZmdnZz+4k2XdmvvP9OowxBhEREREf4md3ASIiIiJlpQAjIiIiPkcBRkRERHyOAoyIiIj4HAUYERER8TkKMCIiIuJzFGBERETE5wTYXcCZ4vF42LlzJ+Hh4TgcDrvLERERkZNgjCEzM5OaNWvi53fs4yyVNsDs3LmTxMREu8sQERGRU7Bjxw5q1659zOWVNsCEh4cD1gcQERFhczUiIiJyMjIyMkhMTPR+jx9LpQ0wJaeNIiIiFGBERER8zImaf6gRr4iIiPgcBRgRERHxOQowIiIi4nMqbRuYk+V2uyksLLS7DJGzyul0HvfyRBGRiq7KBhhjDKmpqaSlpdldishZ5+fnR1JSEk6n0+5SREROSZUNMCXhpXr16oSEhKizO6kySjp53LVrF3Xq1NHPvoj4pCoZYNxutze8xMbG2l2OyFkXFxfHzp07KSoqIjAw0O5yRETKrEqeBC9p8xISEmJzJSL2KDl15Ha7ba5EROTUVMkAU0KHzqWq0s++iPi6Kh1gRERExDcpwIiIiIjPUYAROYFHH32U8847z+4yRETkMAowPuSmm27C4XAwbty4UvM//fTTs96mweFweG8RERG0b9+e//3vf2e1huOxI3R8/PHHdOvWjcjISMLCwmjZsiWPP/44Bw4cOKt1iIiUmTFQmAc5ByBrj3XL3A2ZqZCxCzJ2QvqfkP6HdUvbAWnboSDbtpKr5GXUviwoKIhnnnmGf/zjH0RHR9tay5QpU+jVqxcZGRm88sorXHXVVaxcuZIWLVrYVpMxxpYra/75z3/yzDPPMGLECJ5++mlq1qzJxo0bmTx5Mu+88w733HPPWa9JRKowdyHs3wS7f4Xdv8CBLVCYAwU51n1hLhRmW/cl8zBlf50r/wstrir38k+GjsBgfenlFBTZcjOmbD8wPXr0ICEhgbFjxx53vUWLFnHhhRcSHBxMYmIiw4YNIzvbSsqTJk3i3HPP9a5bcgRn8uTJpV7n4YcfPu5rREVFkZCQQOPGjXniiScoKipi/vz53uU7duzgmmuuISoqipiYGPr168fWrVu9y2+66Sb69+/PY489RlxcHBEREdxxxx0UFBR418nPz2fYsGFUr16doKAgLrjgApYvX+5dvmDBAhwOB7NmzaJt27a4XC7effddHnvsMVavXu09SvTmm28CkJaWxq233up9vYsuuojVq1eXel/jxo0jPj6e8PBwBg8eTF5e3nE/h2XLlvH000/z3HPPMWHCBM4//3zq1avHJZdcwscff8ygQYMA2Lx5M/369SM+Pp6wsDDat2/P119/XWpbr7zyCo0aNSIoKIj4+HiuuurQHwaPx8PYsWNJSkoiODiYVq1a8dFHH3mXHzx4kIEDBxIXF0dwcDCNGjViypQpx61dRHycMdaRkk3zYPFE+OQfMPkCeLomvNIJPh4Mi56HtZ/Cxq9g2yLYuRL2rrOOoGTvtYLMEeHFYd0cfuDwB78A8AsEfyf4u6xbQJC13CY6AgPkFrppPnqOLa+99vEUQpwnvxv8/f15+umnueGGGxg2bBi1a9c+Yp3NmzfTq1cvnnzySd544w327t3L0KFDGTp0KFOmTKFr164MGzaMvXv3EhcXx8KFC6lWrRoLFizgjjvuoLCwkCVLlvDQQw+dVE1FRUX897//BQ71L1JYWEhKSgrJycl89913BAQE8OSTT9KrVy9+/vln73rz5s0jKCiIBQsWsHXrVm6++WZiY2N56qmnABg5ciQff/wxb731FnXr1mX8+PGkpKSwadMmYmJivDU89NBDPPvss9SvX5+goCDuu+8+Zs+e7Q0IkZGRAFx99dUEBwcza9YsIiMjee2117j44ovZsGEDMTExfPDBBzz66KO8/PLLXHDBBbzzzjtMnDiR+vXrH/P9v/fee4SFhXHXXXcddXlUVBQAWVlZXHrppTz11FO4XC7efvtt+vbty/r166lTpw4//vgjw4YN45133uH888/nwIEDfPfdd97tjB07lnfffZfJkyfTqFEjvv32W/72t78RFxdH165deeSRR1i7di2zZs2iWrVqbNq0idzc3JPahyLiI9xFsGs1bFkIWxdZ0zn7jr6uMxzim0P8OVCtMbjCITCk+BYMzlDr/vB5gSHg7xvRwGHKegjAR2RkZBAZGUl6ejoRERGlluXl5bFlyxaSkpIICgoip6DIJwLMTTfdRFpaGp9++inJyck0b96c//73v3z66acMGDDAezTn1ltvxd/fn9dee8373EWLFtG1a1eys7NxuVzExcUxefJkrrrqKlq3bs21117Liy++yK5du1i8eDHdu3cnLS3tmJ39ORwOgoKC8Pf3Jzc3F4/HQ7169VixYgUxMTG8++67PPnkk6xbt87bPqegoICoqCg+/fRTevbsyU033cTnn3/Ojh07vK8zefJkHnjgAdLT08nNzSU6Opo333yTG264AbCCUb169Rg+fDgPPPAACxYsoHv37nz66af069fPW9+jjz7Kp59+yqpVq0p9Bn369GHPnj24XC7v/IYNGzJy5Ehuv/12zj//fFq3bs3LL7/sXd6pUyfy8vJKbetwl156KX/++ecRR3JOxrnnnssdd9zB0KFD+eSTT7j55pv5448/CA8PL7Vefn4+MTExfP311yQnJ3vn33rrreTk5DB16lQuv/xyqlWrxhtvvHHC1/3r74CIVFAeD+xZC1u+tW7bFkN+xl9WckBsAyuoxJ9bfH8ORNYBHxy09Xjf34fzjZh1hgUH+rP28RTbXvtUPPPMM1x00UXcf//9RyxbvXo1P//8M++99553njEGj8fDli1baNasGV26dGHBggX06NGDtWvXctdddzF+/Hh+++03Fi5cSPv27U/YU/Hzzz9Pjx49+P333xkxYgQTJ070HhVZvXo1mzZtOuKLOC8vj82bN3sft2rVqtTrJCcnk5WVxY4dO0hPT6ewsJDOnTt7lwcGBtKhQwfWrVtXarvt2rU74We2evVqsrKyjhg+Ijc311vTunXruOOOO0otT05OLnVq7K9O9n+ArKwsHn30UWbOnMmuXbsoKioiNzeX7du3A3DJJZdQt25d6tevT69evejVqxcDBgwgJCSETZs2kZOTwyWXXFJqmwUFBbRu3RqAO++8kyuvvJKVK1fSs2dP+vfvz/nnn39StYlIBWEM7N9sHWHZ8i1s/Q5y9pdeJygS6l4ASV2gdnuo3gycVa9neQUYrKMJZTmNUxF06dKFlJQURo0axU033VRqWVZWFv/4xz8YNmzYEc+rU6cOAN26deP111/nu+++o3Xr1kRERHhDzcKFC+natesJa0hISKBhw4Y0bNiQKVOmcOmll7J27VqqV69OVlYWbdu2LRWiSsTFxZ3amz6O0NDQE66TlZVFjRo1WLBgwRHLSk7znIrGjRuzaNEiCgsLjzuu0P3338/cuXN59tlnadiwIcHBwVx11VXeNj/h4eGsXLmSBQsW8NVXXzF69GgeffRRli9fTlZWFgAzZ86kVq1apbZbcjSpd+/ebNu2jS+//JK5c+dy8cUXM2TIEJ599tlTfm8ichZk7YHfF8Lv8+H3BZDxZ+nlgaFQN9kKLEldIKEl+J3aP7+ViW99a0sp48aN47zzzqNJkyal5rdp04a1a9fSsGHDYz63a9euDB8+nA8//JBu3boBVqj5+uuvWbx4Mffdd1+ZaunQoQNt27blqaee4sUXX6RNmzZMnz6d6tWrH/cQ4OrVq8nNzSU4OBiAH374gbCwMBITE6lWrRpOp5PFixdTt25dwDqFtHz5coYPH37cepxO5xFXI7Vp04bU1FQCAgKoV6/eUZ/XrFkzli5dyo033uid98MPPxz3tW644QYmTpzIK6+8ctSrjdLS0oiKimLx4sXcdNNNDBgwALAC1eGNmgECAgLo0aMHPXr0YMyYMURFRfHNN99wySWX4HK52L59+3HDZVxcHIMGDWLQoEFceOGFPPDAAwowIhVNQTZsW3IosOz+pfRyfyckdjwUWGq2gQCnLaVWZAowPqxFixYMHDiQiRMnlpr/4IMP0qlTJ4YOHcqtt95KaGgoa9euZe7cuUyaNAmAli1bEh0dzdSpU/niiy8AK8Dcf//9OByOUqdtTtbw4cMZMGAAI0eOZODAgUyYMIF+/frx+OOPU7t2bbZt28Ynn3zCyJEjvY2PCwoKGDx4MA8//DBbt25lzJgxDB06FD8/P0JDQ7nzzjt54IEHiImJoU6dOowfP56cnBwGDx583Frq1avHli1bWLVqFbVr1yY8PJwePXqQnJxM//79GT9+PI0bN2bnzp3MnDmTAQMG0K5dO+655x5uuukm2rVrR+fOnXnvvff49ddfj9uIt2PHjowcOZL77ruPP//8kwEDBlCzZk02bdrE5MmTueCCC7jnnnto1KgRn3zyCX379sXhcPDII4/g8Xi82/niiy/4/fff6dKlC9HR0Xz55Zd4PB6aNGlCeHg4999/PyNGjMDj8XDBBReQnp7O4sWLiYiIYNCgQYwePZq2bdtyzjnnkJ+fzxdffEGzZs3KvB9FpJx5PLBrFWyeZx1p2bEU3AWl10loCfW7QYPukNipSp4SKjNTSaWnpxvApKenH7EsNzfXrF271uTm5tpQ2akbNGiQ6devX6l5W7ZsMU6n0/x1Vy5btsxccsklJiwszISGhpqWLVuap556qtQ6/fr1MwEBASYzM9MYY4zb7TbR0dGmU6dOJ6wFMDNmzCg1z+PxmKZNm5o777zTGGPMrl27zI033miqVatmXC6XqV+/vrntttu8+6Tk/YwePdrExsaasLAwc9ttt5m8vDzvNnNzc83dd9/t3Ubnzp3NsmXLvMvnz59vAHPw4MFSteTl5Zkrr7zSREVFGcBMmTLFGGNMRkaGufvuu03NmjVNYGCgSUxMNAMHDjTbt2/3Pvepp54y1apVM2FhYWbQoEFm5MiRplWrVif8TKZPn266dOliwsPDvZ/5448/7q1ty5Ytpnv37iY4ONgkJiaaSZMmma5du5p77rnHGGPMd999Z7p27Wqio6NNcHCwadmypZk+fXqpz/eFF14wTZo0MYGBgSYuLs6kpKSYhQsXGmOMeeKJJ0yzZs1McHCwiYmJMf369TO///77UWv11d8BEZ9RmG/Mxq+N+XyEMc82MWZMROnbv8415n9DjVnzkTFZe+2utkI53vf34XQVkq7AsM3hV1XJ2aXfAZEzIC8DNs2F32bCxrmlrxZyhh06wlK/O8TUB40Kf1S6CklERORMy9gF67+0br8vBE/hoWVh8dCkNzS9zGrLEuA69nakzBRgRERETpYxVudxG+bAhtlWr7aHq9YYmlxqhZZabX2yHxZfoQAjtinp3l9EpEIryLH6Zdkw2woumbtKL6/dHpr2gSZ9IK6xPTVWQQowIiIif5X+56HAsmUhFB02JlpgqNWWpXEvaNQTwuPtq7MKU4AREREpzIVt31t9s2xeALvXlF4emWgFlia9rF5wA9X43W4KMCIiUvV4PFYHcpu/sULLtiXgzj9sBQckdoDGKdC4t9Vdv64aqlAUYEREpGpI/7P4CEtxD7h/HcU5vCY0uKj4UuduEFrNjirlJCnAiIhI5ZX+B/w6A375GHb+VHpZYCgkXWj1y9Kgu3UFkY6y+AwFGClXDoeDGTNm0L9//wqxHRGpgjJ3w9r/WaFlx2FjmTn8rHGFSjqTq91eYwz5MF2g7mNSU1O5++67qV+/Pi6Xi8TERPr27cu8efPsLu2UPProo5x33nlHzN+1axe9e/c+46//008/cfXVVxMfH09QUBCNGjXitttuY8OGDWf8tUWkHOUcgB+nwFt94V9NYdYDxeHFAXU7Q5/n4L4NcNs8uOhhqNdZ4cXH6QiMD9m6dSudO3cmKiqKCRMm0KJFCwoLC5kzZw5Dhgzht99+s7vEcpOQkHDGX+OLL77gyiuvJCUlhffee48GDRqwZ88ePvzwQx555BGmT59+xmsQkdNQkAPrPoM1H1ltWzxFh5bVbg/nXAHn9IeImraVKGfQWRmZyQaVcTDH3r17m1q1apmsrKwjlh0+YCBgfvrpp1LLADN//nxjzKEBEGfPnm3OO+88ExQUZLp37252795tvvzyS9O0aVMTHh5urr/+epOdne3dTt26dc3zzz9f6nVbtWplxowZ433MXwZ5HDlypGnUqJEJDg42SUlJ5uGHHzYFBQXGGGOmTJligFK3kkEXD99OcnKyGTlyZKnX3bNnjwkICPAOZJiXl2fuu+8+U7NmTRMSEmI6dOjgfb9Hk52dbapVq2b69+9/1OUln2dRUZG55ZZbTL169UxQUJBp3LixeeGFF0qtO3/+fNO+fXsTEhJiIiMjzfnnn2+2bt3qXf7pp5+a1q1bG5fLZZKSksyjjz5qCgsLjTHWAI1jxowxiYmJxul0mho1api77777mHWXF1/9HRAxHo8xO3405rNhxjxVq/QAia9eYMx3zxtzYOsJNyMV18kO5qgjMGB1DV2YY89rB4acVKOxAwcOMHv2bJ566ilCQ0OPWB4VFVXml3700UeZNGkSISEhXHPNNVxzzTW4XC6mTp1KVlYWAwYM4KWXXuLBBx8s87ZLhIeH8+abb1KzZk3WrFnDbbfdRnh4OCNHjuTaa6/ll19+Yfbs2Xz99dcAREZGHrGNgQMHMn78eMaNG4ej+LOaPn06NWvW5MILLwRg6NChrF27lmnTplGzZk1mzJhBr169WLNmDY0aNTpim3PmzGHfvn2MHDnyqHWXfJ4ej4fatWvz4YcfEhsby/fff8/tt99OjRo1uOaaaygqKqJ///7cdtttvP/++xQUFLBs2TJvnd999x033ngjEydO5MILL2Tz5s3cfvvtAIwZM4aPP/6Y559/nmnTpnHOOeeQmprK6tWrT/nzFqm0svfDz9Php3dgz9pD86PrQasb4NwroVpD28qTs08BBqzw8rRNhxj/byc4jwwkf7Vp0yaMMTRt2rTcXvrJJ5+kc+fOAAwePJhRo0axefNm6tevD8BVV13F/PnzTyvAPPzww97pevXqcf/99zNt2jRGjhxJcHAwYWFhBAQEHPeU0TXXXMPw4cNZtGiRN7BMnTqV66+/HofDwfbt25kyZQrbt2+nZk1rP95///3Mnj2bKVOm8PTTTx+xzY0bNwKc8PMMDAzkscce8z5OSkpiyZIlfPDBB1xzzTVkZGSQnp7OZZddRoMGDQBo1qyZd/3HHnuMhx56iEGDBgFQv359nnjiCUaOHMmYMWPYvn07CQkJ9OjRg8DAQOrUqUOHDh2OW5NIleFxW5c8//Q2/PbloYESA4Kg2eXQ5u9Wp3Iab6hKUoDxEcaYct9my5YtvdPx8fGEhIR4w0vJvGXLlp3Wa0yfPp2JEyeyefNmsrKyKCoqOu7w6EcTFxdHz549ee+997jwwgvZsmULS5Ys4bXXXgNgzZo1uN1uGjcuPQZJfn4+sbGxR91mWT7Pl19+mTfeeIPt27eTm5tLQUGBt+FxTEwMN910EykpKVxyySX06NGDa665hho1agCwevVqFi9ezFNPPeXdntvtJi8vj5ycHK6++mpeeOEF6tevT69evbj00kvp27cvAQH61ZQqLPcgLHkFVk2FjD8Oza9xnhVazr0KgqPsqk4qCP2VBOs0zv/ttO+1T0KjRo1wOBwnbKjrV/yfyOFf0IWFhUddNzAw0DvtcDhKPS6Z5/F4Sm37r1/8x9o2wJIlSxg4cCCPPfYYKSkpREZGMm3aNJ577rnjvoejGThwIMOGDeOll15i6tSptGjRghYtWgCQlZWFv78/K1aswN/fv9TzwsLCjrq9krDz22+/kZycfMzXnTZtGvfffz/PPfccycnJhIeHM2HCBJYuXepdZ8qUKQwbNozZs2czffp0Hn74YebOnUunTp3Iysriscce44orrjhi20FBQSQmJrJ+/Xq+/vpr5s6dy1133cWECRNYuHDhEftDpErITIW3+8He4r91wdHQ8lpo/TdIaGFvbVKhKMCA1QblJE7j2CkmJoaUlBRefvllhg0bdkQ7mLS0NKKiooiLiwOsy5Bbt24NwKpVq8qlhri4OHbtOjQKa0ZGBlu2bDnm+t9//z1169bln//8p3fetm3bSq3jdDpxu90nfO1+/fpx++23M3v2bKZOncqNN97oXda6dWvcbjd79uzxnmI6kZ49e1KtWjXGjx/PjBkzjlhe8nkuXryY888/n7vuusu7bPPmzUes37p1a1q3bs2oUaNITk5m6tSpdOrUiTZt2rB+/XoaNjz2ufng4GD69u1L3759GTJkCE2bNmXNmjW0adPmpN6LSKWRtgPevhwO/G71ipvypDXCs8YdkqNQgPEhL7/8Mp07d6ZDhw48/vjjtGzZkqKiIubOncurr77KunXrCA4OplOnTowbN46kpCT27NlTqh3K6bjooot488036du3L1FRUYwePfqIIx6Ha9SoEdu3b2fatGm0b9+emTNnHhEW6tWrx5YtW1i1ahW1a9cmPDwcl8t1xLZCQ0Pp378/jzzyCOvWreP666/3LmvcuDEDBw7kxhtv5LnnnqN169bs3buXefPm0bJlS/r06XPU7f3nP//h6quv5vLLL2fYsGE0bNiQffv28cEHH3jrbtSoEW+//TZz5swhKSmJd955h+XLl5OUlATAli1beP3117n88supWbMm69evZ+PGjd6ANXr0aC677DLq1KnDVVddhZ+fH6tXr+aXX37hySef5M0338TtdtOxY0dCQkJ49913CQ4Opm7duqe0j0R81v7N1pGX9B0QVRcGfWY10BU5ljN/QZQ9KuNl1MYYs3PnTjNkyBBTt25d43Q6Ta1atczll19e6pLhtWvXmuTkZBMcHGzOO+8889VXXx31MuqSS4WNsS5pjoyMLPVaY8aMMa1atfI+Tk9PN9dee62JiIgwiYmJ5s033zzhZdQPPPCAiY2NNWFhYebaa681zz//fKnXycvLM1deeaWJioo65mXUJb788ksDmC5duhzxuRQUFJjRo0ebevXqmcDAQFOjRg0zYMAA8/PPPx/381y+fLm54oorTFxcnHG5XKZhw4bm9ttvNxs3bvTWd9NNN5nIyEgTFRVl7rzzTvPQQw95P5fU1FTTv39/U6NGDeN0Ok3dunXN6NGjjdvt9r7G7Nmzzfnnn2+Cg4NNRESE6dChg3n99deNMcbMmDHDdOzY0URERJjQ0FDTqVMn8/XXXx+35vLgy78DUgntXmfMhMbWpdAT2xiT9ofdFYmNTvYyaocxZ6B1aAWQkZFBZGQk6enpRzQazcvLY8uWLSQlJREUpEOTUvXod0AqjF2r4Z0BkLMfqp8DN34KYdXtrkpsdLzv78PpFJKIiNjjjx/h3SsgLx1qtoa/fQIhMXZXJT5CAUZERM6+rYtg6rVQkAWJnWDgBxB0ZEeWIseiACMiImfXpq9h2kAoyoOkrnD9+xX+SlCpeNR9oYiInD2/zYT3r7fCS6MUuOEDhRc5JVX6CEwlbb8sckL62RdbrPkIPrkdjBua94Mr/gMBTrurEh9VJY/AlPRwmpNj0wCOIjYrKCgAOG4/PiLlauPcQ+Gl5XVw5RsKL3JaquQRGH9/f6KiotizZw8AISEh3tGDRSo7j8fD3r17CQkJ0ZhLcnbsWg0fDLLCS6vrod8rGoBRTluV/etVMvpxSYgRqUr8/PyoU6eOgruceel/WFcbFWZbDXb7TlR4kXJRZQOMw+GgRo0aVK9e/bgDEopURk6n0zvwp8gZk5cO710Nmbsgrhlc+45OG0m5qbIBpoS/v7/aAYiIlDd3oXXaaM9aCIuHgR+qnxcpV/oXTEREypcx8MVw+H0+BIZal0pHJdpdlVQyCjAiIlK+vnsWfnoXHH5w9RSoeZ7dFUklpAAjIiLl5+cP4Jsnrene46Fxir31SKWlACMiIuVj6yL49C5r+vy7ocNt9tYjlZoCjIiInL6962HaDeAptHrZ7fG43RVJJacAIyIipydrD7x3lXXZdO0OMOA19fUiZ5x+wkRE5NQV5Fgd1aVth+gka2TpwGC7q5IqQAFGRERO3RfDYedKCI6Bv30ModXsrkiqCAUYERE5Naveh5+ng8Mfrn0XYhvYXZFUIWUKMG63m0ceeYSkpCSCg4Np0KABTzzxBMYY7zrGGEaPHk2NGjUIDg6mR48ebNy4sdR2Dhw4wMCBA4mIiCAqKorBgweTlZVVap2ff/6ZCy+8kKCgIBITExk/fvxpvE0RESlX+zbBzPus6W6joF5ne+uRKqdMAeaZZ57h1VdfZdKkSaxbt45nnnmG8ePH89JLL3nXGT9+PBMnTmTy5MksXbqU0NBQUlJSyMvL864zcOBAfv31V+bOncsXX3zBt99+y+233+5dnpGRQc+ePalbty4rVqxgwoQJPProo7z++uvl8JZFROS0FOXDRzdbAzTWuxAuvNfuiqQqMmXQp08fc8stt5Sad8UVV5iBAwcaY4zxeDwmISHBTJgwwbs8LS3NuFwu8/777xtjjFm7dq0BzPLly73rzJo1yzgcDvPnn38aY4x55ZVXTHR0tMnPz/eu8+CDD5omTZqcdK3p6ekGMOnp6WV5iyIiciKzRhkzJsKYcfWMSf/T7mqkkjnZ7+8yHYE5//zzmTdvHhs2bABg9erVLFq0iN69ewOwZcsWUlNT6dGjh/c5kZGRdOzYkSVLlgCwZMkSoqKiaNeunXedHj164Ofnx9KlS73rdOnSBafz0KilKSkprF+/noMHDx61tvz8fDIyMkrdRESknG34Cn542Zru/wpE1LS3HqmyyjQa9UMPPURGRgZNmzbF398ft9vNU089xcCBAwFITU0FID4+vtTz4uPjvctSU1OpXr166SICAoiJiSm1TlJS0hHbKFkWHR19RG1jx47lscceK8vbERGRsshMhU/vsKY73gFNettbj1RpZToC88EHH/Dee+8xdepUVq5cyVtvvcWzzz7LW2+9dabqO2mjRo0iPT3de9uxY4fdJYmIVB4eD3xyO+Tsh4QWcIl62hV7lekIzAMPPMBDDz3EddddB0CLFi3Ytm0bY8eOZdCgQSQkJACwe/duatSo4X3e7t27Oe+88wBISEhgz549pbZbVFTEgQMHvM9PSEhg9+7dpdYpeVyyzl+5XC5cLldZ3o6IiJysxc/DloUQGAJXTYEA/b0Ve5XpCExOTg5+f+ke2t/fH4/HA0BSUhIJCQnMmzfPuzwjI4OlS5eSnJwMQHJyMmlpaaxYscK7zjfffIPH46Fjx47edb799lsKCwu968ydO5cmTZoc9fSRiIicQTuWwzdPWdOXToBqjeytR4QyBpi+ffvy1FNPMXPmTLZu3cqMGTP417/+xYABAwBwOBwMHz6cJ598ks8++4w1a9Zw4403UrNmTfr37w9As2bN6NWrF7fddhvLli1j8eLFDB06lOuuu46aNa3GYDfccANOp5PBgwfz66+/Mn36dF588UXuvVeX6omInFW5afDxLWDccO5VcN5AuysSsZTl0qaMjAxzzz33mDp16pigoCBTv359889//rPU5c4ej8c88sgjJj4+3rhcLnPxxReb9evXl9rO/v37zfXXX2/CwsJMRESEufnmm01mZmapdVavXm0uuOAC43K5TK1atcy4cePKUqouoxYROV0ejzEfDLIumX6+hTG5aXZXJFXAyX5/O4w5rBvdSiQjI4PIyEjS09OJiIiwuxwREd+z8m347G7wC4BbvoLabe2uSKqAk/3+1lhIIiJypL3r4cuR1vRFjyi8SIWjACMiIqUVZMNHt0BRLtTvDucPs7sikSMowIiIyCHGwKd3we5fIDQOBrwGfvqqkIpHP5UiInLIt8/C2k/BLxCueQfC40/4FBE7KMCIiIjlt5kw/0lrus9zUDfZ3npEjkMBRkREYPdaa6gAgA63Q9tB9tYjcgIKMCIiVV3OAXj/OijIgqQukPK03RWJnJACjIhIVeYuhA9uhLRtEF0Prn4L/APtrkrkhBRgRESqsjn/B1u/A2cYXPc+hMTYXZHISVGAERGpqla8Cctet6aveB3im9tajkhZKMCIiFRF276Hmfdb090fhqZ97K1HpIwUYEREqpq07TD97+AphOb9ocv9dlckUmYKMCIiVUlBNky7AXL2QUIL6P8KOBx2VyVSZgowIiJVRckwAalrIKSa1WjXGWp3VSKnRAFGRKSqWPT8oWECrn0XohLtrkjklCnAiIhUBduWwDfFwwRcOkHDBIjPU4AREanssvfDR7eAcUOLa6DtTXZXJHLaFGBERCozjwc+vRMyd0JsQ7jsX2q0K5WCAoyISGW2ZBJsnAP+Lrj6TXCF212RSLlQgBERqax2LId5j1nTvcdZl02LVBIKMCIilVHOAfjoZvAUwTlXQNub7a5IpFwpwIiIVDbGwP+GQPoOiKkPfV9UuxepdBRgREQqmx9ehfVfgr/TavcSFGF3RSLlTgFGRKQy+XMFzB1tTac8DTVa2VuPyBmiACMiUlnkpsGHNxUP0tgP2t9qd0UiZ4wCjIhIZWAMfDbUGmk6qi5c/pLavUilpgAjIlIZLPs3rPvcGufo6jchKNLuikTOKAUYERFft3MVfPVPa7rnE1Crja3liJwNCjAiIr6spN2LuwCa9IGOd9hdkchZoQAjIuKrPB6YcQcc3AKRdaDfJLV7kSpDAUZExFct+hdsmGWNc3Tt2xASY3dFImeNAoyIiC/aNA++edKa7vMc1Gxtbz0iZ5kCjIiIr0nbDh/fChhoMwja/N3uikTOOgUYERFfUpgHH9wIuQesoy69x9tdkYgtFGBERHzJrJGw8ycIjoFr3obAILsrErGFAoyIiK9Y+TasfAtwwFX/hag6dlckYhsFGBERX/DnSph5vzV90cPQ4CJ76xGxmQKMiEhFl73favfizocml8IF99pdkYjtFGBERCoyjxs+uRXSd0BMfej/KvjpT7eIfgtERCqyBWNh8zcQEAzXvgvBUXZXJFIhKMCIiFRU62fBtxOs6csnQvw59tYjUoEowIiIVEQHfodP/mFNd/gHtLzG3npEKhgFGBGRimjOw5CfDrU7QM8n7a5GpMJRgBERqWj+WAHrZ4LDD/q9DAFOuysSqXAUYEREKpr5xUdcWl4HcY3trUWkglKAERGpSLYutq468guAbg/aXY1IhaUAIyJSURgD3xQffWlzI0TXs7UckYpMAUZEpKLYPA+2fw8BQdDlAburEanQFGBERCqCw4++tL8VImraW49IBacAIyJSEfw2E3b+BIGh0Hm43dWIVHgKMCIidvO4Yf5T1nSnOyEszt56RHyAAoyIiN1++QT2rIWgSDj/brurEfEJCjAiInZyF8GCp63p8+/WYI0iJ0kBRkTETqunWuMehVSDjnfaXY2Iz1CAERGxS1E+LBxvTV94L7jC7K1HxIcowIiI2GXFW5C+A8JrQLtb7K5GxKcowIiI2KEgB7571pru8gAEBttbj4iPUYAREbHDstchazdE1YHWf7e7GhGfowAjInK25WXA4hes6W6jIMBpazkivkgBRkTkbPvhFcg9CNUaQ8tr7a5GxCeVOcD8+eef/O1vfyM2Npbg4GBatGjBjz/+6F1ujGH06NHUqFGD4OBgevTowcaNG0tt48CBAwwcOJCIiAiioqIYPHgwWVlZpdb5+eefufDCCwkKCiIxMZHx48ef4lsUEalAcg7A95Os6W6jwM/f3npEfFSZAszBgwfp3LkzgYGBzJo1i7Vr1/Lcc88RHR3tXWf8+PFMnDiRyZMns3TpUkJDQ0lJSSEvL8+7zsCBA/n111+ZO3cuX3zxBd9++y233367d3lGRgY9e/akbt26rFixggkTJvDoo4/y+uuvl8NbFhGx0eIXoSAT4ltA8/52VyPiu0wZPPjgg+aCCy445nKPx2MSEhLMhAkTvPPS0tKMy+Uy77//vjHGmLVr1xrALF++3LvOrFmzjMPhMH/++acxxphXXnnFREdHm/z8/FKv3aRJk5OuNT093QAmPT39pJ8jInJG7d1gzBPVjRkTYcxvs+yuRqRCOtnv7zIdgfnss89o164dV199NdWrV6d169b8+9//9i7fsmULqamp9OjRwzsvMjKSjh07smTJEgCWLFlCVFQU7dq1867To0cP/Pz8WLp0qXedLl264HQeatiWkpLC+vXrOXjw4FFry8/PJyMjo9RNRKTCcBfCJ7dDUR7U7w6NU+yuSMSnlSnA/P7777z66qs0atSIOXPmcOeddzJs2DDeeustAFJTUwGIj48v9bz4+HjvstTUVKpXr15qeUBAADExMaXWOdo2Dn+Nvxo7diyRkZHeW2JiYlnemojImfXdv2DnSmvAxn4vg8Nhd0UiPq1MAcbj8dCmTRuefvppWrduze23385tt93G5MmTz1R9J23UqFGkp6d7bzt27LC7JBERy58rYeEz1vSlz0FkLXvrEakEyhRgatSoQfPmzUvNa9asGdu3bwcgISEBgN27d5daZ/fu3d5lCQkJ7Nmzp9TyoqIiDhw4UGqdo23j8Nf4K5fLRURERKmbiIjtCnNhxj/AuOGcAdDiKrsrEqkUyhRgOnfuzPr160vN27BhA3Xr1gUgKSmJhIQE5s2b512ekZHB0qVLSU5OBiA5OZm0tDRWrFjhXeebb77B4/HQsWNH7zrffvsthYWF3nXmzp1LkyZNSl3xJCJS4X39GOzbAGEJ0OdfOnUkUk7KFGBGjBjBDz/8wNNPP82mTZuYOnUqr7/+OkOGDAHA4XAwfPhwnnzyST777DPWrFnDjTfeSM2aNenfvz9gHbHp1asXt912G8uWLWPx4sUMHTqU6667jpo1awJwww034HQ6GTx4ML/++ivTp0/nxRdf5N577y3fdy8icib9vgCWvmpN95sEITG2liNSqZT18qbPP//cnHvuucblcpmmTZua119/vdRyj8djHnnkERMfH29cLpe5+OKLzfr160uts3//fnP99debsLAwExERYW6++WaTmZlZap3Vq1ebCy64wLhcLlOrVi0zbty4MtWpy6hFxFY5B415rpl1yfTnw+2uRsRnnOz3t8MYY+wOUWdCRkYGkZGRpKenqz2MiJx9n/wDfp4G0UlwxyJwhdldkYhPONnvb42FJCJS3tb+zwovDj+44nWFF5EzQAFGRKQ8Ze6Gz4db0xeMgMQOtpYjUlkpwIiIlBdj4LO7IfcAJLSArg/ZXZFIpaUAIyJSXla+BRvngL8TBrwOAc4TP0dETokCjIhIeTiwBWb/nzV90SMQ3/z464vIaVGAERE5XR43zLgDCrOhbmdIHmJ3RSKVngKMiMjpWvIy7PgBnGHQ/1Xw87e7IpFKTwFGROR07N8M85+ypnuNhei69tYjUkUowIiInCqPx7rqqCgP6neH1n+3uyKRKkMBRkTkVK14A7YthsBQ6PuiBmoUOYsUYERETkXaDpg7xpruMUanjkTOMgUYEZGyMga+GA4FWZDYCdrfZndFIlWOAoyISFmtngabvgZ/F1z+EvjpT6nI2abfOhGRssjcDbOLhwjo9iDENba3HpEqSgFGRKQsvrwf8tIgoSWcP8zuakSqLAUYEZGTtfZ/sO4z8AuAfi+Df6DdFYlUWQowIiInI+cAzLzfmu48HGq0tLUckapOAUZE5GTM+Sdk74FqTaDrSLurEanyFGBERE5k49eweirggH6TIMBld0UiVZ4CjIjI8eRnWn2+AHS6ExI72FqOiFgUYEREjufrRyF9B0TXg4setrsaESmmACMicixbF8Py/1jTfSeCM9TeekTESwFGRORoCvOskaYB2gyC+l3trUdESlGAERE5mkXPw4HNEJYAPZ+wuxoR+QsFGBGRv9q3CRb9y5ruPQ6CIu2tR0SOoAAjInI4Y2DmveAugIY9oHl/uysSkaNQgBEROdyaj2DLQggIgkufBYfD7opE5CgUYERESuQehDmjrOkuD0BMkr31iMgxKcCIiJSY9zhk77WGC9BI0yIVmgKMiAjAjuXw4xRr+rJ/QYDT3npE5LgUYERE3EXwxQjAQKsboN4FdlckIiegACMisuw12L0GgqLU54uIj1CAEZGqLf0P+OYpa/qSxyG0mr31iMhJUYARkapt9kNQmA2JnaD13+2uRkROkgKMiFRd62fDus/BL8BquOunP4kivkK/rSJSNRVkw5cPWNPJQyD+HHvrEZEyUYARkapp4XhI3w6RidD1QburEZEyUoARkapn91pYMsmavnQCOEPtrUdEykwBRkSqFo/HGqzRUwRNL4Mmve2uSEROgQKMiFQtq96D7UsgMBR6jbO7GhE5RQowIlJ15B6Erx+1prs9BFGJtpYjIqdOAUZEqo75YyFnnzVYY6c77a5GRE6DAoyIVA2pv8Dyf1vTl44H/0B76xGR06IAIyKVnzFWny/GA837Q/1udlckIqdJAUZEKr81H8H27yEwBHo+aXc1IlIOFGBEpHLLy4CvHramL7xPDXdFKgkFGBGp3L4dD1mpEFMfzr/b7mpEpJwowIhI5bV3PfzwqjXd6xkIcNlbj4iUGwUYEamcjIFZI60edxv3hsY97a5IRMqRAoyIVE7rPoPfF4C/C3o9bXc1IlLOFGBEpPIpyIHZ/2dNd77Hav8iIpWKAoyIVD6L/gUZf0BkHbhghN3ViMgZoAAjIpXL/s2w+EVrutfT4Ayxtx4ROSMUYESkcpk9CtwF0OAiaHqZ3dWIyBmiACMilcf62bBxDvgFQu/x4HDYXZGInCEKMCJSORTmwewHrenku6BaI3vrEZEzSgFGRCqH71+Cg1shvAZ0ecDuakTkDFOAERHft38zfPesNd3zSXCF21uPiJxxCjAi4ts8HvjsbijKg/rd4Nwr7a5IRM4CBRgR8W0rpsC2xRAYAn1fVMNdkSritALMuHHjcDgcDB8+3DsvLy+PIUOGEBsbS1hYGFdeeSW7d+8u9bzt27fTp08fQkJCqF69Og888ABFRUWl1lmwYAFt2rTB5XLRsGFD3nzzzdMpVUQqo/Q/YO4Ya/riMRBdz9ZyROTsOeUAs3z5cl577TVatmxZav6IESP4/PPP+fDDD1m4cCE7d+7kiiuu8C53u9306dOHgoICvv/+e9566y3efPNNRo8e7V1ny5Yt9OnTh+7du7Nq1SqGDx/Orbfeypw5c061XBGpbIyBz4dDQSbU7gAdbrO7IhE5m8wpyMzMNI0aNTJz5841Xbt2Nffcc48xxpi0tDQTGBhoPvzwQ++669atM4BZsmSJMcaYL7/80vj5+ZnU1FTvOq+++qqJiIgw+fn5xhhjRo4cac4555xSr3nttdealJSUk64xPT3dACY9Pf1U3qKIVHSrphkzJsKYx6sZs+c3u6sRkXJyst/fp3QEZsiQIfTp04cePXqUmr9ixQoKCwtLzW/atCl16tRhyZIlACxZsoQWLVoQHx/vXSclJYWMjAx+/fVX7zp/3XZKSop3G0eTn59PRkZGqZuIVFJZew/1+dL1QYhrYm89InLWBZT1CdOmTWPlypUsX778iGWpqak4nU6ioqJKzY+Pjyc1NdW7zuHhpWR5ybLjrZORkUFubi7BwcFHvPbYsWN57LHHyvp2RMQXzXoAcg9CQgtrtGkRqXLKdARmx44d3HPPPbz33nsEBQWdqZpOyahRo0hPT/feduzYYXdJInImrPsCfp0BDn+4fBL4B9pdkYjYoEwBZsWKFezZs4c2bdoQEBBAQEAACxcuZOLEiQQEBBAfH09BQQFpaWmlnrd7924SEhIASEhIOOKqpJLHJ1onIiLiqEdfAFwuFxEREaVuIlLJ5KbBzPus6c7DoOZ5dlYjIjYqU4C5+OKLWbNmDatWrfLe2rVrx8CBA73TgYGBzJs3z/uc9evXs337dpKTkwFITk5mzZo17Nmzx7vO3LlziYiIoHnz5t51Dt9GyTol2xCRKuqrhyErFWIbQteH7K5GRGxUpjYw4eHhnHvuuaXmhYaGEhsb650/ePBg7r33XmJiYoiIiODuu+8mOTmZTp06AdCzZ0+aN2/O3//+d8aPH09qaioPP/wwQ4YMweVyAXDHHXcwadIkRo4cyS233MI333zDBx98wMyZM8vjPYuIL9o8H356B3BYp44CK9ZpbBE5u8rciPdEnn/+efz8/LjyyivJz88nJSWFV155xbvc39+fL774gjvvvJPk5GRCQ0MZNGgQjz/+uHedpKQkZs6cyYgRI3jxxRepXbs2//nPf0hJSSnvckXEFxRkw+fDrOkOt0FdHY0VqeocxhhjdxFnQkZGBpGRkaSnp6s9jIivm/UQLH0VIhPhrh/AFWZ3RSJyhpzs97fGQhKRim37Ulg62Zru+6LCi4gACjAiUpEV5sFnQwED5w2EhhfbXZGIVBAKMCJScX33LOzbAKHVoeeTdlcjIhWIAoyIVEypa2DR89Z0n2chJMbeekSkQlGAEZGKx10E/xsKniJodjk072d3RSJSwSjAiEjFs2QS7FoFQZFw6bN2VyMiFZACjIhULPs2wYKx1nTKWAiPP/76IlIlKcCISMXh8Vgd1hXlQYOL4Lwb7K5IRCooBRgRqThWTIFtiyEwFC57ARwOuysSkQpKAUZEKob0P2DuGGu6xxiIrmtvPSJSoSnAiIj9jIEv7oWCTEjsCO1vtbsiEangFGBExH5rPoSNc8DfCZe/BH7+dlckIhWcAoyI2CtrL8x60JruOhLimthbj4j4BAUYEbHX7Ach9wDEnwudh9tdjYj4CAUYEbHPb1/CLx+Dwx/6TQL/QLsrEhEfoQAjIvbIS4eZ91rT598NNVvbW4+I+BQFGBGxx9zRkLkLYhpAt4fsrkZEfIwCjIicfVu+hRVvWtOXvwSBwbaWIyK+RwFGRM6u7P3w6RBrut1gqNfZ3npExCcpwIjI2eMuhA9uhPTtEF0Pejxqd0Ui4qMUYETk7Jk1ErYtAmc4XD8dgiLsrkhEfJQCjIicHcv/Az++ATjgyv9A9aZ2VyQiPkwBRkTOvC3fwpcjrekeY6BJL3vrERGfpwAjImfWgd+tdi/GDS2vVW+7IlIuFGBE5MzJy4D3r4fcg1CrLfSdCA6H3VWJSCWgACMiZ4bHDZ/cDnt/g/AacO17EBhkd1UiUkkowIjImfHNk7BhFvi74Lr3IKKG3RWJSCWiACMi5e/nD2HRv6zpfi9bp49ERMqRAoyIlK8/V8BnQ63pC0ZAy6vtrUdEKiUFGBEpPxm7YNpAKMqDxr3gokfsrkhEKikFGBEpH4V5MH2gNcJ0XDO44t/g5293VSJSSSnAiMjpMwY+H2adPgqOhuunapgAETmjFGBE5PR9/xL8PB0c/nD1WxBT3+6KRKSSU4ARkdOzcS7MHW1N9xoH9bvaW4+IVAkKMCJy6vZthI8GAwbaDIIOt9ldkYhUEQowInJqctPg/esgPx3qJMOlz2qYABE5axRgRKTsPG74eDDs3wQRteGadyDAaXdVIlKFKMCISNl9PQY2fQ0BwXD9+xAWZ3dFIlLFKMCISNmset+66ghgwKtQo6W99YhIlaQAIyIn748f4fN7rOkuD8A5A+ytR0SqLAUYETk5JcMEuPOhSR/o9n92VyQiVZgCjIicWGEuTLsBslKLhwl4Dfz050NE7KO/QCJyfMZYp412riweJuB9cIXbXZWIVHEKMCJyfEcME5Bkd0UiIgowInIcW76zLpkGDRMgIhWKAoyIHF1mKnx0CxgPtLpewwSISIWiACMiR3IXWWMcZe+B6udAn39pmAARqVAUYETkSPOfhG2LwBkO17wNzhC7KxIRKUUBRkRKWz8LFj1vTfd7Cao1tLceEZGjUIARkUMOboUZ/7CmO96hnnZFpMJSgBERS2EefHAj5KVDrXZwyRN2VyQickwKMCJimTMKdq2G4Bi4+k0IcNpdkYjIMSnAiAj8/AH8+AbggCv+DVGJdlckInJcCjAiVd2e3w6NMN11JDTqYW89IiInQQFGpCrLz4IP/g6FOVC/G3R90O6KREROigKMSFVlDHw+DPZtgPAacMV/wM/f7qpERE6KAoxIVbX8P/DLx8WDNL4JYXF2VyQictIUYESqoj9XwOxR1vQlj0OdTvbWIyJSRgowIlVN9n6YfiN4CqHpZZA8xO6KRETKTAFGpCrxuOHjwZDxB8Q0gP6vaJBGEfFJZQowY8eOpX379oSHh1O9enX69+/P+vXrS62Tl5fHkCFDiI2NJSwsjCuvvJLdu3eXWmf79u306dOHkJAQqlevzgMPPEBRUVGpdRYsWECbNm1wuVw0bNiQN99889TeoYgcsmAs/D4fAkPg2nchKNLuikRETkmZAszChQsZMmQIP/zwA3PnzqWwsJCePXuSnZ3tXWfEiBF8/vnnfPjhhyxcuJCdO3dyxRVXeJe73W769OlDQUEB33//PW+99RZvvvkmo0eP9q6zZcsW+vTpQ/fu3Vm1ahXDhw/n1ltvZc6cOeXwlkWqqPWz4dsJ1nTfiRDf3N56REROg8MYY071yXv37qV69eosXLiQLl26kJ6eTlxcHFOnTuWqq64C4LfffqNZs2YsWbKETp06MWvWLC677DJ27txJfHw8AJMnT+bBBx9k7969OJ1OHnzwQWbOnMkvv/zifa3rrruOtLQ0Zs+efVK1ZWRkEBkZSXp6OhEREaf6FkUqhwNb4PWu1jhHHW6HSyfYXZGIyFGd7Pf3abWBSU9PByAmJgaAFStWUFhYSI8eh3rybNq0KXXq1GHJkiUALFmyhBYtWnjDC0BKSgoZGRn8+uuv3nUO30bJOiXbOJr8/HwyMjJK3UQEKMy1OqvLS4fa7aHnU3ZXJCJy2k45wHg8HoYPH07nzp0599xzAUhNTcXpdBIVFVVq3fj4eFJTU73rHB5eSpaXLDveOhkZGeTm5h61nrFjxxIZGem9JSZqLBcRjIGZ90HqGgipBle/pUEaRaRSOOUAM2TIEH755RemTZtWnvWcslGjRpGenu697dixw+6SROy34k1Y9R44/OCqNyCylt0ViYiUi4BTedLQoUP54osv+Pbbb6ldu7Z3fkJCAgUFBaSlpZU6CrN7924SEhK86yxbtqzU9kquUjp8nb9eubR7924iIiIIDg4+ak0ulwuXy3Uqb0ekcvpzBcwaaU1f9AjU72pvPSIi5ahMR2CMMQwdOpQZM2bwzTffkJSUVGp527ZtCQwMZN68ed5569evZ/v27SQnJwOQnJzMmjVr2LNnj3eduXPnEhERQfPmzb3rHL6NknVKtiEiJ5C9Hz4YBO4CaNIHLhhhd0UiIuWqTEdghgwZwtSpU/nf//5HeHi4t81KZGQkwcHBREZGMnjwYO69915iYmKIiIjg7rvvJjk5mU6drK7Ke/bsSfPmzfn73//O+PHjSU1N5eGHH2bIkCHeIyh33HEHkyZNYuTIkdxyyy188803fPDBB8ycObOc375IJeRxwye3QvoOiKkPA15VZ3UiUvmYMgCOepsyZYp3ndzcXHPXXXeZ6OhoExISYgYMGGB27dpVajtbt241vXv3NsHBwaZatWrmvvvuM4WFhaXWmT9/vjnvvPOM0+k09evXL/UaJyM9Pd0AJj09vUzPE/F58540ZkyEMU/EG7Nrjd3ViIiUycl+f59WPzAVmfqBkSpp/Wx4/1presDr0Opae+sRESmjs9IPjIhUEMbA8v/C9L9Zj9vfqvAiIpXaKV2FJCIVSGEufHEvrJ5qPW7WF1KetrcmEZEzTAFGxJcd2GL1spu6xurr5eIx0PkeNdoVkUpPAUbEV22cCx/fCnlpVi+7V72hvl5EpMpQgBHxNR4PfDseFowDDNRqC9e8DZG1T/hUEZHKQgFGxJfkHoRPboeNX1mP2w2GXmMhQL1Qi0jVogAj4it2/Wy1dzm4FQKC4LLn4bwb7K5KRMQWCjAivmDV+/DFcCjKg6i6cO27UKOl3VWJiNhGAUakIsvPhFkPWiNKAzTqCVe8DsHR9tYlImIzBRiRiuqPFfDxYDi4xbpEuuuD0GUk+Kn/SRERBRiRisbjhkX/gvljwbghMtE66lL3fLsrExGpMBRgRCqStB0w4x+wbbH1+Nwroc+/IDjK1rJERCoaBRiRiuKXT6yGunnp4AyDS5+FVtepV10RkaNQgBGx218b6tZqB1f+G2Lq21uXiEgFpgAjYqe/NtS98D6rsa5/oN2ViYhUaAowInYoyodFL1hDAniK1FBXRKSMFGBEzrZt38Pnw2HfeuvxOVdYveqqoa6IyElTgBE5W3IOwNdjYOXb1uPQOOg1zrrSSA11RUTKRAFG5EwzBtZ8CLNHQc4+a16bQXDJY+pRV0TkFCnAiJxJB36HL+6F3+dbj+OawmUvQN1kW8sSEfF1CjAiZ0JRAXw/Eb6dYA3A6O+Crg/A+fdAgNPu6kREfJ4CjEh52/6D1Uh37zrrcVJXq5FubANbyxIRqUwUYETKS2YqfP0YrJ5qPQ6JhZSx0PIaNdIVESlnCjAip6soH3541TpdVJBlzWv9N7jkCQiJsbc2EZFKSgFG5FQZAxvmwJxRVmNdgFptofd4qN3O3tpERCo5BRiRU7FvI8x+CDZ9bT0Oi4cej0LL68DPz9bSRESqAgUYkbLIS4eF42HpZGsIAL9ASL4LLrwfgiLsrk5EpMpQgBE5GR6PNVr0vMcge681r3EvSHlaVxeJiNhAAUbkePIzYfU0WPY67NtgzYttZA0B0KiHvbWJiFRhCjAiR7NvoxVaVr0PBZnWPFcEdB0JHf6hzuhERGymACNSwuO2ripa9vqhrv8BqjWGDrdDy2vVzkVEpIJQgBHJOQA/vQPL/wNp24tnOqBJbyu41O+mjuhERCoYBRipunavhR9ehjUfWeMVAQRFQZsbof2tEF3X1vJEROTYFGCkajEGNs+DJS/D5m8OzY9vAR1vh3OvAmeIffWJiMhJUYCRqqEoH37+wAouJYMsOvyg6WXQ6S6o00mniUREfIgCjFRu2fvhx//Csn9D9h5rXmCodZqo4z8gJsne+kRE5JQowEjltG+jdbRl9fuH2rdE1LJCS5tBEBxla3kiInJ6FGCkcvlzBXz7HKyfeWhejVaQfDec0x/8A20rTUREyo8CjPg+Y2DbYvj22dL9tzTuDecPhbqd1b5FRKSSUYAR32WMNRr0t8/Cjh+seQ5/aHkNXDAC4prYW5+IiJwxCjDiezwe+O0L+O5Z2LXamufvhNZ/g873QHQ9W8sTEZEzTwFGfIe7CH75GBb9C/b+Zs0LDIF2t0DyUIioYW99IiJy1ijASMVmDPy5Etb9D36ZAenFXf27Iq2O5zreCaGx9tYocgbkFrhZuf0gbetGExTob3c5IhWOAoxUPB4P7FgKa/8H6z6HjD8OLQupBsl3WV39B0XaV6PIGZKZV8g7P2zjv99tYX92AfViQxh3ZUs61VdQFzmcAoxUDO4i2LYI1n5mtW/J2n1omTMMGvWE5pdDoxR19S+V0sHsAqYs3sKb328lI68IAD8HbN2fw3Wv/8DAjnV4qHdTwoPUFYAIKMCInXIPwtZFsGEO/DYTcg8cWuaKtEaDbn45NLgIAoPtq1N8WkGRhx0Hc9i6L5st+7LZuj+bPw7m4grwIzI4kIigQOs+uOQ+oNT8yJBAXAFn7hTOnow8/v3d77y3dDs5BW4AGlYP465uDejepDoTvlrP1KXbeW/pdr75bQ9PD2hB96bVz1g9Ir7CYYwxdhdxJmRkZBAZGUl6ejoRERF2lyMABTmwfQlsWQi/Lyy+guiwH7+QWGjaB5r1g6QuEOC0rVSpGIrcHtbuyiCv0IPDASW9+Vjd+ji88xwOBw4gLbfQG1QODytuz+n9mQt3BRAb5iQm1ElsmItqYU5iQ13Fj51UCzs0HRPiJMDf74Tb3HEgh9e+3cwHP/5BQZEHgHNqRjC0e0NSzknAz+9Q30Xfb97HqE/WsG1/DgD9z6vJ6L7nEBOq3xGpfE72+1sBRs4cd6HVM+7vC63QsmMZeApLr1OtCdTvag2qWLcz+OugYFVXUORh8aZ9zPplF3PX7uZgTuGJn3QCIU5/6sWGklQtlHrVQkiMDqHQY8jILSQjt5D03EIy8orvc4u8jzNyCzmV7BMVEkhscdiJLQ42saEuYsOcRIc4WbhhL5/+9CdFxRtvWzeaod0b0q1JHI5jdLqYW+DmX3PX899FW/AYiA118ujl53BZyxrHfI6IL1KAUYCxR34WrP/Sutx5y3dQmF16eURtK7AkdbWOsujSZwHyCt0s3LCX2b+k8vXa3WTmF3mXRQZbYQCs43XGmOJ7MBjrvvivWKjrUFCxwop1Xz3cdUpf8h6PITOviH3Z+RzILmB/Vj77sgoOTWcXcCCrgP3Z+ezPKuBgTkGZAs8FDasxpHtDOtWPOen6Vu1I48GPfmb97kwALmkez5P9zyU+Iui4z8srdJOZV0ReoZtqYS6Cnad3WiynoIit+3L4fV8WezLyaRwfTqvESLXRkdOmAKMAc/a4C2HzfFjzgdWWpTDn0LKQWKh34aHQElNf3fpXANv2ZzNt+Q7ScgopdHsocnso9BiK3B6K3MY7Xej2UOg2BPg5ik+hlJw+cRIT5qJa8VGGmFAn0SGBJ3XqpER2fhHz1+9h1ppU5q/f423/AVA93EXKOQn0bpFAh3oxZdqundweQ1qOFXD2FQcb73RWSQgqoHqEi8EXJNG6TvQpvU5BkYeX52/ilQWbKHQbwoMCuLJNbfKL3GTkFVlHlvKKyMyzjihl5BV6T1OViAl1UisqmJpRQdSKCim+D6ZWdDA1o4KJDXXi9hj+OJjLln3ZbN6b5T01t2VfNrvS846oy+GAxtXDaV0nitZ1omhTJ5oGcWGlToeJnIgCjALMmWUM/LEcfv4Afp0BOfsOLYtOsrrzb3oZxJ8Lfr7x5VMVbN+fw0vfbOSTn/487XYhf+VwWEdLnP5++Dms9inHuncAfxzMJf+wL9VaUcH0OjeB3ucm0KZOtL70TsJvqRk8+NHPrP4j/aTWdzjA6e9X6nM/FleAH26P8Z7mOprokECSqoUSG+Zi3a4M/jiYe8Q64a4AzqsTRes60bSuE8V5taOIVtsdOQ4FGAWYM2PvButIy5oP4eDWQ/ND4+CcK6zgUqutjrJUMDsO5DDpm018vPIP7xdS18ZxtK0bTYC/A6e/HwF+DgL8/Qj0dxDg50dggB+BxfMK3R72/+V0yaF769TJqfwlqRcbQq9za9D73ARa1o5UW45T4PYYPvxxB5v3ZhERFEh4UAARwYGEBwUS4Z227sOcATgckJFXxJ8Hc9mZlsvO9Fz+PJjLn2nWbWdaLnsy8737MyjQj3qxodSPKzk1F0ZStVDqVws9Iojsyczjp+1p/LQ9jZXbD/LzH2nkFR4ZlmpGBnFOrUjOqRnBuTUjOadWBAkRQdr/AijAKMCUF3cR/LEMNsyGDV/B3nWHlgWGQrPLoMU1UL+bGuBWQH+m5TLpm018+OMOb3Dp0jiO4T0a0eYUT18cjdtjOFh86qTIbfAYq22Kp7i9ivXY4DFWuxKDdQqjUfUwfWlVQAVFHlLT8/D3d1AjIuiUj4YVuT38lprJTzvS+GnbQX7akcaWfdlHXTc21EnzmhGcWxxsWtaKok6s+nyqihRgFGBOXc4Ba5TnDXOs+7y0Q8v8AqDBxdaRlia9wRlqW5lybDvTcnl5/iY++HEHhW7rV/yChtUYcUkj2taNsbk6qcoy8gpZuzODX3dm8OvOdH79M4NNe7OOekrzkubx/PPSZtSrpr8zVYkCjALMyTMG9qw9dJTlj2VgDjvsGxwNDS+BxinQ8GLrsVQIxhhyC92k5xaSlmPdvlyzi+nLd1Dgtvbh+Q1iGXFJY9rXU3CRiimv0M1vqZlWoNmZwa9/pvPLzgzcHkOgv4ObOycx9KKGROgKpypBAUYB5uiy91kjOe9db932rYc960p33Q9W49tGPaFxL6jdDvw0mNzZkFvgZl9WfvGt5LJda9oKKdb94beSIyx/1al+DCN6NKajxtARH7RxdyZPzFzHtxv2AtYppvt6NuHa9on4q4F3paYAU9UDTNZe2L3mUFApCSs5+4++fkCQ1Y6lUU/rFpV4VsutLPKL3Gzfn0N6biHZBW5yC4rIKXAX36zp3OLH2QVFVh8jWVZj2H1Z+aUuJS6LAD+H1e19cCB1Y0O4rUt9zm9QrZzfncjZZYxhwfq9PDFzLb/vtdrONE0IZ3Tf5vr5rsQUYKpKgDEGMlNh1yqra/5dq2HnKsjceeznRNWFuCbWrVrxfUKLSjPeUF6hm/3ZBbgC/HAF+BEU6E9gOfcjkpFXyKY9WWzak8Xmkvu9WWw/kHNKPbcezhngR1yYy9tFfbUwq6+V6JBAb0iJDHZ6x+mJCg4kxOmvxrBSaRW6PbyzZBsvfL3BO9Blyjnx/N+lzagbq/YxlY0CTGUMMMZA+h+HwsrO4vvsPUdZ2WF1Gle9WXFYaWrdxzaqFKM5G2PYl1XA5r1WcPh9b7Z3+o+DuUdc0uvv5/CGmaAAP1yB/t7HzgA/nMWXDwf6W5cPl3rs74czwI/cArcVWvZmsTcz/5i1lYybE+wMIMTpT4jTn+BAf0JdAQQ7/QkJLJ7nDCAsKIC44qBSMsZOmCtAYUTkKA5mF/D81xt4b+l23B6D09+Pmy+ox9861iUiKJBgp/X7LL5NAaYyBJj8LNj5k9Vh3J8rrPu/tlUBcPhZAaVGq+LbeZBwLrjCy60UYwwHcwrZlZ5LRm5xD5/FPX5m5lk9fZb0+pmZb80LCvQnLtxF9XAXceEu4sJcVI8IIi7MehwT6jziXHZeods7Bk16bkmPooXeRqrbD+RYQWVPlvc/saMJ9Hccs21IeUmICKJB9VAaxoXRsHoYDaqH0TAujLhT7LZeRE7Oht2ZPPHFWr7buO+IZQF+juJ/HIr/gXD5ExJo/fMQ6vInMthJVPGRy+gQp/coZlRI8fy/jD5ujKHA7aGgyOqVuqDImi5wu8kvnuf2eHB7oMjjwVNyb4y3S4EijylukOzn/SfG+w+N06ovxFX+R4p9VaUIMC+//DITJkwgNTWVVq1a8dJLL9GhQ4eTeq7PBRiPx2qj8sePhwLLnrWlrwYCPI4AMiIakhl9Llkx55ATew75Mc3xc4V4OyAr6Zgs0N8PV6AfQQH+BBUfcThefw75RW7+OJjL9gM57Ci+bT+Qw/YDuew4kENW/rEDw6nw93MQG+okLCjACkG5hSfVQ2gJhwMSo0NoEBdK/bgwGsSF0SAulAbVw4gNdWIMFLg95Bd6yCtyk1do/cHJK3STV+ghv8i6LyzuMt/6Y+ShsOQPlfvQspLu9OuXhJW4UI35ImIjYwzz1+9h/Oz1bN6bVa7/sAQH+uPv5ygOKif/N+l0Bfo7CA70J9jpj7/DgZ+fAz+HA38/B37FPVn7F8/z8wN/h4OwoACiQqxR0KNDncSEBFr3odbAoTHF00GBvnMhhs8HmOnTp3PjjTcyefJkOnbsyAsvvMCHH37I+vXrqV69+gmfX+ECTH4mZOz03tzpf1J08A886TtxZO4kIGMbAYVZRzwtlWqscDfgJ09DfvI05BeTRD6n3g23syTUHHYKxRXgx4HsAlIz8k7Ym2q1MBeRwVavnqV7/QwgIiiweH4AYa4Acgrc7MnMZ2/xbU9mHnszratq9mcfu+dWh4PibQUQWfw6EUFW249a0cFWUKkeSr3YUJ/6pRSRM6fQ7TlmY3nvdH6R1d1AbgHpOYWk5RZyMOfQdNpJDMYZ4OewTjt7TzVb/zT6F4eLI26HzT9U46E6cwrc5T6sx9GEOv1JjAmhbmwI9WJDqRsbSt1Y63GNyOAKdWWXzweYjh070r59eyZNmgSAx+MhMTGRu+++m4ceeuiEzz9TASZ1x2Yy9++kMCcdd04G7rwMPHkZmPxMHPmZ+BVk4VeYiX9hNoFFWYQWHiC6aC8hJueE2842Ln72NGCVacCq4sCyB6vPlQA/B/ERQSREBuEK8KOo+AhBkad48L3iowQlg/KVHEnIL/IcdyyTvyr5Ia9TfEs87L52dHC5BYai4q7p92bmk5VfRHhQcVgp7u5c4+CIyNnm8RiyCopIyy7EbYw3pDiLLwgI9Pcr9y/6klNUuYcFm7xCD26PdfrJYwxuj9Wbtcdj9WbtLpnvNmTlF3Egu8AaRDSngIPZhRwoHt6j5P5ER6ec/n7UjgmmbkwIdWNDiQl1EhzoT1DxP7tBgf7Fj/0JdpaedyaO7pzs93eF7Pu9oKCAFStWMGrUKO88Pz8/evTowZIlS476nPz8fPLzDzWszMjIOCO17Zw+nDZZ357SczNMCLtMDKklN2LY64glMzCOTFcN8qMaUD0qlITIIDpHBnFlcWBJiAyiWqjrtLrzLjl1crRTKPlFbqJCnCRGBxMT6jwr7TcC/P2IjwgiPiLojL+WiMjJ8PNzeI/4ni0OhwNXgD+uAH+izsD1FcZYIWdPZj7b9+ewdX822/bnsG1/NtuKmwoUuD38vje7+FL1vWXa/vPXtmJA69rlX/hJqJABZt++fbjdbuLj40vNj4+P57fffjvqc8aOHctjjz12xmsrDKrG7qxYcv1CyPMLId8vlIKAUAoDQikKCMMdGIbHFQ7OcBzOMExoHETWIjCyJsHhkYS5AqjvCqBVUAChroCz0mgrwN+PAH8/Ql0VcneLiMgZ4nA4CA+yBvdsEBd2xHK3x7ArPbc41OSw7UA2GbmF5BVaR4VK2g/mFnrIL3STW1j8uMBNXpGHoAD7TuNXmm+0UaNGce+993ofZ2RkkJhY/p2xdRw6pdy3KSIiYgd/Pwe1o0OoHR1C54Z2V1M2FTLAVKtWDX9/f3bvLn3J8O7du0lISDjqc1wuFy6X62yUJyIiIjarkBedO51O2rZty7x587zzPB4P8+bNIzk52cbKREREpCKokEdgAO69914GDRpEu3bt6NChAy+88ALZ2dncfPPNdpcmIiIiNquwAebaa69l7969jB49mtTUVM477zxmz559RMNeERERqXoqbD8wp6vCdWQnIiIiJ3Sy398Vsg2MiIiIyPEowIiIiIjPUYARERERn6MAIyIiIj5HAUZERER8jgKMiIiI+BwFGBEREfE5CjAiIiLicxRgRERExOdU2KEETldJB8MZGRk2VyIiIiInq+R7+0QDBVTaAJOZmQlAYmKizZWIiIhIWWVmZhIZGXnM5ZV2LCSPx8POnTsJDw/H4XCU23YzMjJITExkx44dGmOpgtI+qvi0jyo+7aOKr7LuI2MMmZmZ1KxZEz+/Y7d0qbRHYPz8/Khdu/YZ235ERESl+oGpjLSPKj7to4pP+6jiq4z76HhHXkqoEa+IiIj4HAUYERER8TkKMGXkcrkYM2YMLpfL7lLkGLSPKj7to4pP+6jiq+r7qNI24hUREZHKS0dgRERExOcowIiIiIjPUYARERERn6MAIyIiIj5HAaaMXn75ZerVq0dQUBAdO3Zk2bJldpdUZX377bf07duXmjVr4nA4+PTTT0stN8YwevRoatSoQXBwMD169GDjxo32FFsFjR07lvbt2xMeHk716tXp378/69evL7VOXl4eQ4YMITY2lrCwMK688kp2795tU8VVz6uvvkrLli29HaElJycza9Ys73Ltn4pn3LhxOBwOhg8f7p1XVfeTAkwZTJ8+nXvvvZcxY8awcuVKWrVqRUpKCnv27LG7tCopOzubVq1a8fLLLx91+fjx45k4cSKTJ09m6dKlhIaGkpKSQl5e3lmutGpauHAhQ4YM4YcffmDu3LkUFhbSs2dPsrOzveuMGDGCzz//nA8//JCFCxeyc+dOrrjiChurrlpq167NuHHjWLFiBT/++CMXXXQR/fr149dffwW0fyqa5cuX89prr9GyZctS86vsfjJy0jp06GCGDBnifex2u03NmjXN2LFjbaxKjDEGMDNmzPA+9ng8JiEhwUyYMME7Ly0tzbhcLvP+++/bUKHs2bPHAGbhwoXGGGt/BAYGmg8//NC7zrp16wxglixZYleZVV50dLT5z3/+o/1TwWRmZppGjRqZuXPnmq5du5p77rnHGFO1f490BOYkFRQUsGLFCnr06OGd5+fnR48ePViyZImNlcnRbNmyhdTU1FL7KzIyko4dO2p/2SQ9PR2AmJgYAFasWEFhYWGpfdS0aVPq1KmjfWQDt9vNtGnTyM7OJjk5WfunghkyZAh9+vQptT+gav8eVdrBHMvbvn37cLvdxMfHl5ofHx/Pb7/9ZlNVciypqakAR91fJcvk7PF4PAwfPpzOnTtz7rnnAtY+cjqdREVFlVpX++jsWrNmDcnJyeTl5REWFsaMGTNo3rw5q1at0v6pIKZNm8bKlStZvnz5Ecuq8u+RAoyInHFDhgzhl19+YdGiRXaXIn/RpEkTVq1aRXp6Oh999BGDBg1i4cKFdpclxXbs2ME999zD3LlzCQoKsrucCkWnkE5StWrV8Pf3P6Jl9+7du0lISLCpKjmWkn2i/WW/oUOH8sUXXzB//nxq167tnZ+QkEBBQQFpaWml1tc+OrucTicNGzakbdu2jB07llatWvHiiy9q/1QQK1asYM+ePbRp04aAgAACAgJYuHAhEydOJCAggPj4+Cq7nxRgTpLT6aRt27bMmzfPO8/j8TBv3jySk5NtrEyOJikpiYSEhFL7KyMjg6VLl2p/nSXGGIYOHcqMGTP45ptvSEpKKrW8bdu2BAYGltpH69evZ/v27dpHNvJ4POTn52v/VBAXX3wxa9asYdWqVd5bu3btGDhwoHe6qu4nnUIqg3vvvZdBgwbRrl07OnTowAsvvEB2djY333yz3aVVSVlZWWzatMn7eMuWLaxatYqYmBjq1KnD8OHDefLJJ2nUqBFJSUk88sgj1KxZk/79+9tXdBUyZMgQpk6dyv/+9z/Cw8O95+MjIyMJDg4mMjKSwYMHc++99xITE0NERAR33303ycnJdOrUyebqq4ZRo0bRu3dv6tSpQ2ZmJlOnTmXBggXMmTNH+6eCCA8P97YbKxEaGkpsbKx3fpXdT3ZfBuVrXnrpJVOnTh3jdDpNhw4dzA8//GB3SVXW/PnzDXDEbdCgQcYY61LqRx55xMTHxxuXy2Uuvvhis379enuLrkKOtm8AM2XKFO86ubm55q677jLR0dEmJCTEDBgwwOzatcu+oquYW265xdStW9c4nU4TFxdnLr74YvPVV195l2v/VEyHX0ZtTNXdTw5jjLEpO4mIiIicErWBEREREZ+jACMiIiI+RwFGREREfI4CjIiIiPgcBRgRERHxOQowIiIi4nMUYERERMTnKMCIiIiIz1GAEZFKqVu3bgwfPtzuMkTkDFGAEZEzavLkyYSHh1NUVOSdl5WVRWBgIN26dSu17oIFC3A4HGzevPksVykivkYBRkTOqO7du5OVlcWPP/7onffdd9+RkJDA0qVLycvL886fP38+derUoUGDBnaUKiI+RAFGRM6oJk2aUKNGDRYsWOCdt2DBAvr160dSUhI//PBDqfndu3fH4/EwduxYkpKSCA4OplWrVnz00UeltvvLL7/Qu3dvwsLCiI+P5+9//zv79u07Zh0zZ84kMjKS9957z/taHTp0IDQ0lKioKDp37sy2bdvK982LyBmjACMiZ1z37t2ZP3++9/H8+fPp1q0bXbt29c7Pzc1l6dKldO/enbFjx/L2228zefJkfv31V0aMGMHf/vY3Fi5cCEBaWhoXXXQRrVu35scff2T27Nns3r2ba6655qivP3XqVK6//nree+89Bg4cSFFREf3796dr1678/PPPLFmyhNtvvx2Hw3HmPwwRKRcBdhcgIpVf9+7dGT58OEVFReTm5vLTTz/RtWtXCgsLmTx5MgBLliwhPz+fbt260bx5c77++muSk5MBqF+/PosWLeK1116ja9euTJo0idatW/P00097X+ONN94gMTGRDRs20LhxY+/8l19+mX/+8598/vnndO3aFYCMjAzS09O57LLLvKermjVrdrY+DhEpBwowInLGdevWjezsbJYvX87Bgwdp3LgxcXFxdO3alZtvvpm8vDwWLFhA/fr1ycrKIicnh0suuaTUNgoKCmjdujUAq1evZv78+YSFhR3xWps3b/YGmI8++og9e/awePFi2rdv710nJiaGm266iZSUFC655BJ69OjBNddcQ40aNc7gpyAi5UkBRkTOuIYNG1K7dm3mz5/PwYMHvUdCatasSWJiIt9//z3z58/noosuIisrC7DarNSqVavUdlwuF2BdxdS3b1+eeeaZI17r8BDSunVrVq5cyRtvvEG7du1KnSKaMmUKw4YNY/bs2UyfPp2HH36YuXPn0qlTp3J//yJS/hRgROSs6N69OwsWLODgwYM88MAD3vldunRh1qxZLFu2jDvvvJPmzZvjcrnYvn27N+j8VZs2bfj444+pV68eAQHH/jPWoEEDnnvuObp164a/vz+TJk0qtbx169a0bt2aUaNGkZyczNSpUxVgRHyEGvGKyFnRvXt3Fi1axKpVq0oFk65du/Laa69RUFBA9+7dCQ8P5/7772fEiBG89dZbbN68mZUrV/LSSy/x1ltvATBkyBAOHDjA9ddfz/Lly9m8eTNz5szh5ptvxu12l3rdxo0bM3/+fD7++GNvx3Zbtmxh1KhRLFmyhG3btvHVV1+xceNGtYMR8SE6AiMiZ0X37t3Jzc2ladOmxMfHe+d37dqVzMxM7+XWAE888QRxcXGMHTuW33//naioKNq0acP//d//Adapp8WLF/Pggw/Ss2dP8vPzqVu3Lr169cLP78j/y5o0acI333zjPRIzcuRIfvvtN9566y32799PjRo1GDJkCP/4xz/OzochIqfNYYwxdhchIiIiUhY6hSQiIiI+RwFGREREfI4CjIiIiPgcBRgRERHxOQowIiIi4nMUYERERMTnKMCIiIiIz1GAEREREZ+jACMiIiI+RwFGREREfI4CjIiIiPic/we5P1rpz1vhgwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = pd.read_csv('../../data/Ebola/Ebola_in_SL_Data.csv', index_col='Weeks')\n", "data.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recipe\n", "\n", "#### Step 1: Construct an 'error' function\n", "We'll begin by constructing a function which takes the model parameters that we intend to vary, and returns the sum of the squared error between the model's prediction and the reported data.\n", "\n", "Our optimizer will interact with our parameter set through an ordered list of values, so our function will need to take this list and unpack it before we can pass it into our model.\n", "\n", "With `pysd` we can ask directly for the model components that we're interested in, at the timestamps that match our data." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "157977495.47574666" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def error(param_list):\n", " #unpack the parameter list \n", " population, contact_frequency = param_list\n", " #run the model with the new parameters, returning the info we're interested in\n", " result = model.run(params={'total_population':population,\n", " 'contact_frequency':contact_frequency},\n", " return_columns=['population_infected_with_ebola'],\n", " return_timestamps=list(data.index.values))\n", " #return the sum of the squared errors\n", " return sum((result['population_infected_with_ebola'] - data['Cumulative Cases'])**2)\n", "\n", "error([10000, 10])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 2: Suggest a starting point and parameter bounds for the optimizer\n", "The optimizer will want a starting point from which it will vary the parameters to minimize the error. We'll take a guess based upon the data and our intuition.\n", "\n", "As our model is only valid for positive parameter values, we'll want to specify that fact to the optimizer. We know that there must be at least two people for an infection to take place (one person susceptible, and another contageous) and we know that the contact frequency must be a finite, positive value. We can use these, plus some reasonable upper limits, to set the bounds." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "susceptible_population_guess = 9000\n", "contact_frequency_guess = 20\n", "\n", "susceptible_population_bounds = (2, 50000)\n", "contact_frequency_bounds = (0.001, 100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 3: Minimize the error with an optimization function\n", "We pass this function into the optimization function, along with an initial guess as to the parameters that we're optimizing. There are a number of optimization algorithms, each with their own settings, that are available to us through this interface. In this case, we're using the L-BFGS-B algorithm, as it gives us the ability to constrain the values the optimizer will try." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " fun: 22200248.07232056\n", " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([2.60749363e+00, 3.06777631e+03])\n", " message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 93\n", " nit: 10\n", " njev: 31\n", " status: 0\n", " success: True\n", " x: array([8.82124175e+03, 8.20469365e+00])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res = scipy.optimize.minimize(error, [susceptible_population_guess,\n", " contact_frequency_guess],\n", " method='L-BFGS-B',\n", " bounds=[susceptible_population_bounds,\n", " contact_frequency_bounds])\n", "res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Result\n", "If we run the simulation with the parameters suggested by the optimizer, we see that the model follows the general behavior of the data, but is too simple to truly capture the correct shape of the curve." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(2, 9000, 'RMSE: 7.5% of Max')" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHHCAYAAABeLEexAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACRZElEQVR4nOzdd3hTZRvA4V+6d0uBDlYpe++9p5U9VYaCLFHhU4aAqAwHorjAATgBEQQRUQTZU7aA7CmUTQctXXQn7/fHoYHQFpqSEto893XlSnLOm3Oek6QnT991dEophRBCCCGEDbOzdgBCCCGEENYmCZEQQgghbJ4kREIIIYSweZIQCSGEEMLmSUIkhBBCCJsnCZEQQgghbJ4kREIIIYSweZIQCSGEEMLmSUIkhBBCCJsnCZGN0Ol0TJ061ezXXbhwAZ1Ox/z58x9Y9uzZszzxxBN4e3uj0+n4/fffmT9/PjqdjgsXLpi978fF888/T+nSpa0dhjDD1q1b0el0/Prrrxbb5tSpU9HpdBbb3vPPP4+Hh0eOyub27zevmHNeyI2Mz2/r1q33Lfew55eszlmPiqW/T49CQTif348kRI9QxpdJp9OxY8eOTOuVUpQsWRKdTkfnzp2tEOHDGThwIEePHmXatGksXLiQevXqZVlu9uzZOT6RJiYmMnXq1AeeGHMjMjKSV199lUqVKuHq6oqfnx8NGjRgwoQJJCQkWHx/j0rp0qXz5fcnJzJ+RLK7hYWFWTvEx8bixYuZOXOmtcN4rOX0nJVbeXn+ym/yw/fRwdoB2CIXFxcWL15Ms2bNTJZv27aNK1eu4OzsbKXIci8pKYndu3fz5ptvMnLkSOPy5557jj59+pgc0+zZsylSpAjPP//8A7ebmJjI22+/DUCrVq0sFm90dDT16tUjLi6OwYMHU6lSJaKiojhy5Ahz5szhpZdeMv73/u2332IwGCy2b/Hw5syZk2Xtio+Pz6MP5hFISkrCwcG80/XixYs5duwYo0aNsng8QUFBJCUl4ejoaPFtmyOr80tOZXfOsqT7nb/eeustXn/99TzZ7+MoL7+PliIJkRV07NiRZcuW8fnnn5uc5BYvXkzdunW5ceOGFaPLncjISCDzD5K9vT329vZWiOj+vv/+ey5dusTOnTtp0qSJybq4uDicnJyMzy150ldKkZycjKura6Z1ycnJODk5YWcnFbcP0rt3b4oUKWLtMB4ZFxcXa4cAQHp6OgaDAScnp8cipoc5v2R3znpUHBwczE5yRd6SM68V9O3bl6ioKDZs2GBclpqayq+//kq/fv2yfM2tW7cYO3YsJUuWxNnZmYoVK/Lxxx+jlDIpl5KSwujRoylatCienp507dqVK1euZLnNq1evMnjwYPz9/XF2dqZq1ar88MMPZh/P1KlTCQoKAmDcuHHodDpjn5t725xLly7N8ePH2bZtm7GZI7uanwsXLlC0aFEA3n77bWP5u/tSbN68mebNm+Pu7o6Pjw/dunXj5MmTD4z53Llz2Nvb06hRo0zrvLy8TE72WfUhMhgMzJw5k6pVq+Li4oK/vz/Dhw/n5s2bJuUymq/WrVtHvXr1cHV15euvvzb2kViyZAlvvfUWxYsXx83Njbi4OKKjo3nttdeoXr06Hh4eeHl50aFDBw4fPvzA4zLHTz/9RN26dXF1dcXX15c+ffpw+fLlTOWWLVtmLFekSBGeffZZrl69alImoz/M1atX6d69Ox4eHhQtWpTXXnsNvV6fq/fOEvR6PW+88QYBAQG4u7vTtWvXXB9jVubNm0ebNm3w8/PD2dmZKlWqMGfOHLNizMl7du/3Pj4+nlGjRlG6dGmcnZ3x8/Ojffv2HDx4ENBqI1avXs3FixeNfzd3f4cjIiIYMmQI/v7+uLi4ULNmTRYsWGCyz4x+Qh9//DEzZ86kbNmyODs7c+LEiWz7EJ06dYqnn36aokWL4urqSsWKFXnzzTeN6y9evMjLL79MxYoVcXV1pXDhwjz11FO57pOSVZ+WjL+5HTt20KBBA1xcXChTpgw//vijscz9zlmQ83NjcnIyU6dOpUKFCri4uBAYGEjPnj05d+7cA89fWfUhSk9P59133zW+16VLl+aNN94gJSXFpFxOjhEgLS2Nt99+m/Lly+Pi4kLhwoVp1qyZyW9Pdo4fP06bNm1wdXWlRIkSvPfee1nWlP/xxx906tSJYsWK4ezsTNmyZXn33XdNvsP3+z6mpqYyefJk6tati7e3N+7u7jRv3pwtW7Y8MEZLk/TUCkqXLk3jxo35+eef6dChAwBr1qwhNjaWPn368Pnnn5uUV0rRtWtXtmzZwpAhQ6hVqxbr1q1j3LhxXL16lc8++8xYdujQofz000/069ePJk2asHnzZjp16pQphvDwcBo1aoROp2PkyJEULVqUNWvWMGTIEOLi4syq1uzZsyc+Pj6MHj2avn370rFjx2w7i86cOZP//e9/eHh4GE+U/v7+WZYtWrSosfmqR48e9OzZE4AaNWoAsHHjRjp06ECZMmWYOnUqSUlJfPHFFzRt2pSDBw/etyN0UFAQer2ehQsXMnDgwBwfa4bhw4czf/58Bg0axCuvvEJoaChffvkl//77Lzt37jSpVTp9+jR9+/Zl+PDhDBs2jIoVKxrXvfvuuzg5OfHaa6+RkpKCk5MTJ06c4Pfff+epp54iODiY8PBwvv76a1q2bMmJEycoVqyY2fHea9q0aUyaNImnn36aoUOHEhkZyRdffEGLFi34999/jf81Zxxj/fr1mT59OuHh4cyaNYudO3ealAMt+QgJCaFhw4Z8/PHHbNy4kU8++YSyZcvy0ksv5eq9y050dHSmZQ4ODpn+2582bRo6nY4JEyYQERHBzJkzadeuHYcOHTLW0plzjPeaM2cOVatWpWvXrjg4OPDnn3/y8ssvYzAYGDFixAOPI6fv2b1efPFFfv31V0aOHEmVKlWIiopix44dnDx5kjp16vDmm28SGxvLlStXjOeHjL/JpKQkWrVqxX///cfIkSMJDg5m2bJlPP/888TExPDqq6+a7GvevHkkJyfzwgsv4OzsjK+vb5Y/jEeOHKF58+Y4OjrywgsvULp0ac6dO8eff/7JtGnTAPjnn3/YtWsXffr0oUSJEly4cIE5c+bQqlUrTpw4gZub2wPfs5z477//6N27N0OGDGHgwIH88MMPPP/889StW5eqVave95yV03OjXq+nc+fObNq0iT59+vDqq68SHx/Phg0bOHbsGO3atbvv+SsrQ4cOZcGCBfTu3ZuxY8eyd+9epk+fzsmTJ1mxYoVZxwha0jV9+nSGDh1KgwYNiIuLY//+/Rw8eJD27dtnG0dYWBitW7cmPT2d119/HXd3d7755pssa7bnz5+Ph4cHY8aMwcPDg82bNzN58mTi4uL46KOPAO77fYyLi+O7776jb9++DBs2jPj4eL7//ntCQkLYt28ftWrVysEnbiFKPDLz5s1TgPrnn3/Ul19+qTw9PVViYqJSSqmnnnpKtW7dWimlVFBQkOrUqZPxdb///rsC1HvvvWeyvd69eyudTqf+++8/pZRShw4dUoB6+eWXTcr169dPAWrKlCnGZUOGDFGBgYHqxo0bJmX79OmjvL29jXGFhoYqQM2bN+++x5ZR7qOPPsrymENDQ43Lqlatqlq2bHnf7WWIjIzMFHuGWrVqKT8/PxUVFWVcdvjwYWVnZ6cGDBhw3+2GhYWpokWLKkBVqlRJvfjii2rx4sUqJiYmU9mBAweqoKAg4/O///5bAWrRokUm5dauXZtpeVBQkALU2rVrTcpu2bJFAapMmTLG9zpDcnKy0uv1JstCQ0OVs7Ozeuedd+57XBn7vPv7c68LFy4oe3t7NW3aNJPlR48eVQ4ODsblqampys/PT1WrVk0lJSUZy61atUoBavLkycZlAwcOVECm+GrXrq3q1q1rfG7Oe5eVKVOmKCDLW8WKFY3lMt7f4sWLq7i4OOPyX375RQFq1qxZZh9jxr7vdu9np5RSISEhqkyZMvc9DqVy/p4ppTL9DXh7e6sRI0bcd/udOnUy+d5mmDlzpgLUTz/9ZFyWmpqqGjdurDw8PIzvV8bftJeXl4qIiDDZRlbnhRYtWihPT0918eJFk7IGg8H4OKv3a/fu3QpQP/74o3FZxue3ZcuW+x5jVueXjL+57du3G5dFREQoZ2dnNXbs2EzHcO85K6fnxh9++EEB6tNPP80UV8Yx3+/8de/3KeP8PXToUJNyr732mgLU5s2bzT7GmjVr3vdckJ1Ro0YpQO3du9dk+97e3pne76w+0+HDhys3NzeVnJxsXJbd9zE9PV2lpKSYLLt586by9/dXgwcPNjv2hyFNZlby9NNPk5SUxKpVq4iPj2fVqlXZNpf99ddf2Nvb88orr5gsHzt2LEop1qxZYywHZCp3b22PUorly5fTpUsXlFLcuHHDeAsJCSE2NtZY9f64un79OocOHeL555/H19fXuLxGjRq0b9/e+F5kx9/fn8OHD/Piiy9y8+ZN5s6dS79+/fDz8+Pdd9/N1BR5t2XLluHt7U379u1N3ru6devi4eGRqao3ODiYkJCQLLc1cODATP91OTs7G/sR6fV6oqKi8PDwoGLFihb5XH777TcMBgNPP/20SfwBAQGUL1/eGP/+/fuJiIjg5ZdfNmlC7NSpE5UqVWL16tWZtv3iiy+aPG/evDnnz583Pjf3vcvO8uXL2bBhg8lt3rx5mcoNGDAAT09P4/PevXsTGBho/H7k5hjvdvdnFxsby40bN2jZsiXnz58nNjY2R8fyoPcsKz4+Puzdu5dr167laB93++uvvwgICKBv377GZY6OjrzyyiskJCSwbds2k/K9evUyNv1kJzIyku3btzN48GBKlSplsu7uZqG736+0tDSioqIoV64cPj4+Fj3nVKlShebNmxufFy1alIoVKz7wfTXn3Lh8+XKKFCnC//73v0zbyc1w+ozv5JgxY0yWjx07FiDTdzEnx+jj48Px48c5e/as2bE0atSIBg0amGy/f//+mcre/ZnGx8dz48YNmjdvTmJiIqdOnXrgvuzt7Y19Ng0GA9HR0aSnp1OvXr1H/jskTWZWUrRoUdq1a8fixYtJTExEr9fTu3fvLMtevHiRYsWKmZzYASpXrmxcn3FvZ2dH2bJlTcrd3UQD2skrJiaGb775hm+++SbLfUZEROTquB6VjGO+99hAe1/WrVvHrVu3cHd3z3YbgYGBzJkzh9mzZ3P27FnWrVvHhx9+yOTJkwkMDGTo0KFZvu7s2bPExsbi5+eX5fp737vg4OBsY8hqncFgYNasWcyePZvQ0FCTtvjChQtnu62cOnv2LEopypcvn+X6jCar+73HlSpVyjR1hIuLS6YfzkKFCpn0DTL3vctOixYtctSp+t5j1Ol0lCtXztjnxNxjvNfOnTuZMmUKu3fvJjEx0WRdbGws3t7e9319Tt6zrMyYMYOBAwdSsmRJ6tatS8eOHRkwYABlypS57+tAO+by5ctn6rx/7/kkw/2+vxkyfoSrVat233JJSUlMnz6defPmcfXqVZN/PHKaQObEvUkZ5Ox9NefceO7cOSpWrGixjtEZ5+9y5cqZLA8ICMDHxyfT55KTY3znnXfo1q0bFSpUoFq1ajz55JM899xz9222y4ilYcOGmZZn9Xdy/Phx3nrrLTZv3kxcXJzJupx+pgsWLOCTTz7h1KlTpKWlGZfn5LtnSZIQWVG/fv0YNmwYYWFhdOjQ4ZGNdsho/3/22Wez7T/zoD+YgkSn01GhQgUqVKhAp06dKF++PIsWLco2ITIYDPj5+bFo0aIs19/7A5dVu/v91r3//vtMmjSJwYMH8+677+Lr64udnR2jRo2yyPB/g8GATqdjzZo1WY7QyelkgffKyWgfc9+7x9m5c+do27YtlSpV4tNPP6VkyZI4OTnx119/8dlnn+Xos8rtCKmnn36a5s2bs2LFCtavX89HH33Ehx9+yG+//Wbsl2gp9/v+mut///sf8+bNY9SoUTRu3Ng4IWKfPn0sOrVFdu/r/Wp+4fE4N+a0diknx9iiRQvOnTvHH3/8wfr16/nuu+/47LPPmDt3brbnN3PExMTQsmVLvLy8eOeddyhbtiwuLi4cPHiQCRMm5Ogz/emnn3j++efp3r0748aNw8/PD3t7e6ZPn865c+ceOkZzSEJkRT169GD48OHs2bOHpUuXZlsuKCiIjRs3Eh8fb1JLlFEdmTFaIigoCIPBYPzPJcPp06dNtpcxAk2v19OuXTtLHlKOmFOdnF3ZjGO+99hAe1+KFCly39qh7JQpU4ZChQpx/fr1bMuULVuWjRs30rRpU4v+WGT49ddfad26Nd9//73J8piYGIsMNS9btixKKYKDg6lQoUK25e5+j9u0aWOy7vTp08b15u47L9+7e93bVKCU4r///jP+qD3MMf7555+kpKSwcuVKk//WH9XomMDAQF5++WVefvllIiIiqFOnDtOmTTMmRPf72zly5AgGg8Gkluje84k5Mmqmjh07dt9yv/76KwMHDuSTTz4xLktOTiYmJsbsfeYFc86NZcuWZe/evaSlpWU7EMCcc13G+fvs2bPG2jrQOnnHxMTk6nMB8PX1ZdCgQQwaNIiEhARatGjB1KlT75sQBQUFZdnMdu/5duvWrURFRfHbb7/RokUL4/LQ0NBMr83uvfj1118pU6YMv/32m0mZKVOmPPDYLE36EFmRh4cHc+bMYerUqXTp0iXbch07dkSv1/Pll1+aLP/ss8/Q6XTGE2DG/b2j1O6dHdTe3p5evXqxfPnyLE9gGfNz5BV3d/ccnwAzRp3cWz4wMJBatWqxYMECk3XHjh1j/fr1dOzY8b7b3bt3L7du3cq0fN++fURFRWVZNZzh6aefRq/X8+6772Zal56e/tAnd3t7+0z/yS5btixHw8BzomfPntjb2/P2229n2o9SiqioKADq1auHn58fc+fONRn2u2bNGk6ePJnl6MUHyev37l4//vgj8fHxxue//vor169fN/6tPMwxZvyHfm+zT1Z9mSxJr9dnaorw8/OjWLFiJsfg7u6eZZNFx44dCQsLM/knLD09nS+++AIPDw9atmxpdkxFixalRYsW/PDDD1y6dMlk3d3vT1bf7S+++CLTNAPWYs65sVevXty4cSPTeRnuHHN256+sZJyz7j1ff/rppwC5+nvL+FvO4OHhQbly5TIN488qlj179rBv3z7jssjIyEw1u1n9DaSmpjJ79uxM28zu+5jVNvbu3cvu3bvvG2NekBoiK8vJkO8uXbrQunVr3nzzTS5cuEDNmjVZv349f/zxB6NGjTL2GapVqxZ9+/Zl9uzZxMbG0qRJEzZt2sR///2XaZsffPABW7ZsoWHDhgwbNowqVaoQHR3NwYMH2bhxY5bDmi2lbt26zJkzh/fee49y5crh5+eX6b/zDK6urlSpUoWlS5dSoUIFfH19qVatGtWqVeOjjz6iQ4cONG7cmCFDhhiH3Xt7ez/wuk8LFy5k0aJF9OjRg7p16+Lk5MTJkyf54YcfcHFx4Y033sj2tS1btmT48OFMnz6dQ4cO8cQTT+Do6MjZs2dZtmwZs2bNyrY/WE507tyZd955h0GDBtGkSROOHj3KokWLctQ/JMN///3He++9l2l57dq16dSpE++99x4TJ07kwoULdO/eHU9PT0JDQ1mxYgUvvPACr732Go6Ojnz44YcMGjSIli1b0rdvX+OQ9NKlSzN69Gizj81S792vv/6aZdNe+/btTaZx8PX1pVmzZgwaNIjw8HBmzpxJuXLlGDZsGMBDHeMTTzyBk5MTXbp0Yfjw4SQkJPDtt9/i5+d33xrGhxUfH0+JEiXo3bs3NWvWxMPDg40bN/LPP/+Y1LzUrVuXpUuXMmbMGOrXr4+HhwddunThhRde4Ouvv+b555/nwIEDlC5dml9//ZWdO3cyc+bMTH0Vc+rzzz+nWbNm1KlThxdeeIHg4GAuXLjA6tWrOXToEKB9txcuXIi3tzdVqlRh9+7dbNy40SJ94ywlp+fGAQMG8OOPPzJmzBj27dtH8+bNuXXrFhs3buTll1+mW7du9z1/3atmzZoMHDiQb775xtgUtW/fPhYsWED37t1p3bq12cdSpUoVWrVqRd26dfH19WX//v3G6RruZ/z48SxcuJAnn3ySV1991TjsPqN2MUOTJk0oVKgQAwcO5JVXXkGn07Fw4cIsmyaz+z527tyZ3377jR49etCpUydCQ0OZO3cuVapUefSXUHqUQ9ps3d3D7u8nq2HT8fHxavTo0apYsWLK0dFRlS9fXn300UcmQ1qVUiopKUm98sorqnDhwsrd3V116dJFXb58Ocuhn+Hh4WrEiBGqZMmSytHRUQUEBKi2bduqb775xlgmL4bdh4WFqU6dOilPT08FPHAI/q5du1TdunWVk5NTpuPYuHGjatq0qXJ1dVVeXl6qS5cu6sSJE/fdnlJKHTlyRI0bN07VqVNH+fr6KgcHBxUYGKieeuopdfDgQZOy9w67z/DNN9+ounXrKldXV+Xp6amqV6+uxo8fr65du2Ysk90Q+IxhxcuWLcu0Ljk5WY0dO1YFBgYqV1dX1bRpU7V7927VsmXLHE1XkDEkN6vbkCFDjOWWL1+umjVrptzd3ZW7u7uqVKmSGjFihDp9+rTJ9pYuXapq166tnJ2dla+vr+rfv7+6cuVKpvfI3d09UyxZDVXP6XuXlfsNu+euYdoZ7+/PP/+sJk6cqPz8/JSrq6vq1KlTpmHhOT3GrI5l5cqVqkaNGsrFxUWVLl1affjhh8bh2Hd/57Niznt29/c+JSVFjRs3TtWsWVN5enoqd3d3VbNmTTV79myT1yQkJKh+/fopHx8fBZh8h8PDw9WgQYNUkSJFlJOTk6pevXqmv/Hs/qbvXnfva44dO6Z69OihfHx8lIuLi6pYsaKaNGmScf3NmzeN+/Xw8FAhISHq1KlTKigoSA0cONBY7mGH3Wf1N3fv38/9ji8n50altCHnb775pgoODjaW6927tzp37pyxTHbnr6w+57S0NPX2228bt1eyZEk1ceJEk+Hr5hzje++9pxo0aKB8fHyUq6urqlSpkpo2bZpKTU3N9Np7HTlyRLVs2VK5uLio4sWLq3fffVd9//33md7vnTt3qkaNGilXV1dVrFgxNX78eLVu3bpMn19230eDwaDef/99FRQUpJydnVXt2rXVqlWrsj3v5iWdUg/oZSaEEEIIUcBJHyIhhBBC2DxJiIQQQghh8yQhEkIIIYTNk4RICCGEEDZPEiIhhBBC2DxJiIQQQghh82RixhwwGAxcu3YNT0/PXF3FWAghhBCPnlKK+Ph4ihUrlumCxveShCgHrl27RsmSJa0dhhBCCCFy4fLly5QoUeK+ZSQhyoGMqewvX76Ml5eXlaMRQgghRE7ExcVRsmTJHF2SRhKiHMhoJvPy8pKESAghhMhnctLdRTpVCyGEEMLmSUIkhBBCCJsnCZEQQgghbJ4kREIIIYSweZIQCSGEEMLmSUIkhBBCCJsnCZEQQgghbJ4kREIIIYSweZIQCSGEEMLmSUIkhBBCCJsnCZEQQgghbJ4kREIIIYSweXJxVyGEEEJYlkEPaYmQngLKAEoB6vY9dz2+615nD97FrRayJERCCCGEuL/EaAg/rt0iTsCtG1rCk5YEabdu3ydB6u3H+hTz9+ERAK+dtnzsOSQJkRBCCCE0+jS4cfZ28nPsThIUf+0hN6wDne6u+yyWOTg/5D4ejiREQgghhK1KjIYLf0Podri0FyJPgSEt67I+QeBfDfyrgFcxcHQHR1dwcgNHN+2xcdnteweXuxKgx5skREIIIYStSI6DS7u1BCh0G4QdzVzGyRP8q951qwZ+lcHF69HH+whJQiREfnH2LIwYAXv3QlwcrFgB3btbOyohxOMsLQku7bmdAG2Ha/+C0puWKVoZgltA6aYQWAt8SuWbWh1LkmH3Bdn8+dqXOuPm4ADFi8Pzz8PVq5nLt2qllStfPuvtbdhwZ1u//mq67uhR6N0bgoLAxUXbT/v28MUXpuVKlzaN6e7bk0/m7jgz4s7q5uj44Nc//3zWr61UybRcTAz07w+FCkGZMvD995m3tX8/uLlBaGjujuV+Bg7U3udp02DhQqhXL+tyFy7cOYb33su6TP/+2noPD8vHKYSwHoMerh6Evz+BBV3gg1KwsDvs+BSu7teSoULBUGcg9Poexp6BEXug4wyo0g0KBdlkMgRSQ2Qb3nkHgoMhORn27NESpR074NgxLXm5m4sL/Pcf7NsHDRqYrlu0SFufnGy6fNcuaN0aSpWCYcMgIAAuX9b2NWsW/O9/puVr1YKxYzPHWaxY7o7vzTdh6FDTZbduwYsvwhNP5Gwbzs7w3Xemy7y9TZ+/9hps3Qpvv629R8OGQeXK0KSJtl4peOUVGDVKe78tKSkJdu/WjnXkyJy9xsUFfv4Z3nrLdPmtW/DHH5k/eyFE/hQdCue3wPmtWi1Q0k3T9Z7FoEzL27VAzcGnpFXCfNxJQmQLOnS4U5swdCgUKQIffggrV8LTT5uWLVsW0tO1H9K7E6LkZK2JplMnWL7c9DXTpmnJwz//gI+P6bqIiMzxFC8Ozz770Idl1L595mU//aTd9++fs204ODw4plWrYMYMGDBAe37kCPz5552EaNEiuHgR3ngjZ/s0R2Skdn/v+3s/HTvCb7/B4cNQs+ad5X/8AampWo3c5s0WDVMI8Qgkx8K5zVoCdG4LxFw0Xe/spSU+ZVtDmVZQuJzN1vqYQ5rMbFHz5tr9uXNZr+/bF5YuBYPhzrI//4TExMwJVMZ2qlbN+sfazy93MaalwalTcP167l6/eDG4u0O3bjl/jV6v9c3JTlKS1lyWwddXe09Aq3V5/XWYPt38Zqh//9WSVi8v7bVt22q1axmmTtWaIgHGjdNObKVLP3i7jRtrNVWLF5suX7RIS4Z8fTO/5o8/tKS3WDGt1qxsWXj3Xe29yXDyJLi63kkMM+zYAfb2MGFCTo5aCGGOuGuw71tY2ANmlIVlz8OB+VoyZOcAQU2h9ZswZAOMD4W+i6HBMChSXpKhHJIaIlt04YJ2f/eP+9369dN+hLduhTZttGWLF2s/1FklOEFBWnPOsWNQrdqD95+WBjduZF7u7q790ILWx6lyZa3fzPz5D97m3SIjtf5OzzyjbTMnEhO1hCQxUXtf+vbVatHuTm7q14dPP9X6Fp0/D2vXwrffauvef1+r+XruOfNiPX5cS1C9vGD8eK3P09dfa/2itm2Dhg2hZ08t2Rw9WourY8ecJ119+2q1ZR98oJ0Ub9yA9eu1Pkhr12YuP3++tu0xY7T7zZth8mQtUfzoI61M5cpakjRunNZvrGtXLSF8/nntvXnnHfPeAyFEZkppQ+BPrdZu1w6ari9SAcq1gzKtIagJOEt/wIemxAPFxsYqQMXGxlo7FPPMm6cUKLVxo1KRkUpdvqzUr78qVbSoUs7O2vO7tWypVNWq2uN69ZQaMkR7fPOmUk5OSi1YoNSWLdo2ly2787r165Wyt9dujRsrNX68UuvWKZWamjmmoCB1e6L2zLfp0++UCw3Vlg0caP5xf/GF9tq//spZ+ddfV2rCBKWWLlXq55+1fYJSTZsqlZZ2p9yRI0qVKHEn3l69lNLrlTp/XilXV6V27zY/1u7dtff23Lk7y65dU8rTU6kWLe4sy3g/Pvrowdu8u+yxY9rjv//W1n31lVIeHkrduqUdp7u76WsTEzNvb/hwpdzclEpOvrNMr1eqWTOl/P2VunFDqREjlHJwUOqff3J44EKITPTpSl3crdS6N5WaVUupKV533byV+q69UjtmKhV51tqR5hvm/H5LDZEtaNfO9Hnp0lqtQYkS2b+mXz+tFmD2bG1Emb099OgBBw5kLtu+vVZDNH06rFunPZ4xA4oW1Toqd+1qWr5hw6xHP909uq106buueWOmxYu1fWfVtygr06ebPu/TBypU0Dow//qr9hygenVt6PuxY1qNTbly2vKxY6FXL2jUSOuz8/bbWo3KoEEwaVL21dV6vVZb0727NmotQ2Cg9v5/+622Ha+HmPujalWoUUPrE9asmfbedOumjYTLSkYNHUB8PKSkaDVYX3+tNWFm9EWys9Nqk2rW1Jr79u/XOm9nN/JNCJG15DitP9CZdXB2PSTeVXtu76z1AarUCSp2AI9cdkEQOSIJkS346ivtBz42Fn74AbZv1/qH3E+fPtqoqjVrtD4nnTuDp2f25evX15KB1FStE++KFfDZZ1qTyqFDUKXKnbJFimRO0izl/HktIRs5UusonVujR2vJzMaNdxIi0EZm3f2jv3mzltScPq3d+vTRkofSpbXmqpIltcQoK5GRWhNdxYqZ11WurPXhunxZS2oeRr9+8Mkn2jHt2nX/Tt/Hj2uJzebNmftTxcaaPi9bVmtaHTdOayqdNOnh4hTCVkSHwpm12u3CTtOZoV28ocKTWhJUtq00hT1CkhDZggYN7vyId++u1RT066f9gGfXFyUwUOvH8sknsHNn5pFl2XFy0pKj+vW1JGzQIFi2DKZMscSRPFhGB+Kcji7LjqsrFC4M0dHZl9Hr4dVXtc7UxYtrNWpNmtxJgIYP15LJ7BKiR6VvX5g4UZsmoHDh7KciiImBli21Gql33tESHhcXOHhQ6yh9dyf7DOvXa/fXrkFUlDblghDClD4druyD02u0mqAb91zAtHA5LQmq8CSUagT2OZg/TVicJES2xt5eayJq3Rq+/FL7Mc9Ov37aMH0fH60jr7kykrDcjhTLjcWLtR/yRo0ebjvx8VoH5KJFsy8zZ45W7rXXtOfXrpnOpVSsWNYTYGYoWlRrujqdxdWdT53SmqVKWmC+kFKloGlTrZP8Sy9lX3O2dauW1Pz2G7RocWd5dpNMzp2rdV6fNk37Tg0fro1SE0JAzCVtSHzG8PjkmDvr7BygVOM7SVCRctaKUtxFEiJb1KqVVms0c6Y2iWB2E/T17q012VSsqNX8ZGfLljuzRd/tr7+0+6yahB4kLU0bzu/trdVW5cS//2pDwu/XdJMx1UDZstp9crK2r3ubA999V+vDlN3s2dHRWq3X3Ll33j9/f+2yGhlOnrx/jYm9vVZb88cf2si/jKH04eFaYtes2cP1H7rbe+9pn9Mzz9w/HjDtu5WaqvUju1doqNZU1quX1gRXuLA2EeaPP2Yeji+ELUiO0y6SmpEERd8zrYlrISj/BFQI0ZrCXH2sEqbIniREtmrcOHjqKa1j7IsvZl3G21vrI/Ig//uf1hemRw9t2HVqqtZXZelS7Uf+3iajq1fvTJx4Nw+PO9fmys2w+0WLtPv7NZe1bavdZ0w9EBYGtWtrzUoZl+pYt05L5p58Mvt5jCZN0jpZP/XUnWW9emlNTS+9pE1F8PXX2jD9+3nvPa2WpVkzePllrfbm66+1zswzZjzwkHOsZUvtdj9NmmhTDgwcqM24rdNpw/Pv7dyuFAwerDUrzpmjLRs+XGtWffVVrX9YbmcdFyK/0KfD1QPaDNHnNsOV/abXCNPZQ4l62rD4sq2hRH2ws7devOKBJCGyVT17arUkH3+s9S2xf4g/1I8/1voJ/fUXfPONlhCVKqX9wL/1VuYJGw8dynq+nqCg3F+s1GCAJUugTh3zaqR8fLQO4xs2wIIFWr+gcuW0eYVee01rtrrX0aPa6Lm7a4NAS5DmzdOSyPh47fhfeOH++69aFf7+W+vjM326dhwNG2oJY8OGOT8OSyhcWJuNe+xY7XMrVEibvbttWwgJuVPuiy+05rXly02bFL//XutcPWwYrF79aGMX4lEw6OHCDji2HE6uzHyJDN+yWvJTtg2UbqZ1kBb5hk6p3I5tth1xcXF4e3sTGxuLl6WaMIQQQjz+DAa4vBeO/wbHf4dbd12OyLUQBLe8fYmM1tqFUcVjxZzfb6khEkIIIe6mlDYz9LHf4PgKiLtrcIRrIajcFar10mqBpBmswJCESAghhACIPANHlmhNYjcv3Fnu7KXNC1StlzZRogyLL5AkIRJCCGG7UhK0WqB/F2pNYxkc3bTZoav21K4Z5pjNaFxRYEhCJIQQwrYoBVf+gYM/aslQaoK2XGcP5dtDjWe04fFOObw49GNMKYXeoDAoMNx+rFcKZQC9cZ0yrsu4WGPGMoPStqG4vcyg3SsFioz7O2WUunvZ3ctvx4O2Uhnju7MdJwc7GpUpbI23CZCESAghhK1IiITDP8O/P5nOFu1bFmo/C7X6gadlZltP1xtISEknPjmdhJR0ktL0JKfqSU7Xk5Rq0J7fviXdtTxVryc13UCaXpGabiBVbyBNb7i9zHB7mSJNb0BvUKQbDKTrFekGLaFJv7087fZzvSH/jJvy83Rm35t5dFmnHJCESAghRMF2YSfsma1dO8yQri1zdIMq3aHOc9qs0dlchDk13UD0rVRuJKQQdSuVqIQU7XFCKjGJacSnpBGfnH77lmZMgBJT9Vlu73Fkb6fDXqdDpwM7nQ672/c6HdjZ6dCR8Vxbl1FOB+huv2+6LJbrAHSYPs8oi87kLdfpdBR2v88EwI+AJERCCCEKroM/wspXMDbSFK8LtZ+Dar1IcXDnys0kLp+J5HJ0IpdvJnHlZiI34lO5cUtLemKT0u67+QdxdrDD08UBVyd7XB3tcbl90x7bmSxzub3M0d4OJ3s7nBxuP3aww9FeZ7JMu+mwt9PhYGeHvZ3O9Lm9Dge7jOdaMpOR+NjZod3rdNjZZZ0I2iJJiIQQQhRMe7+GNeMBuFy8I9v8BvBvSjEuH0jk8oZ/CItLzjQRe1bs7XT4ujtR2N2Jop7OFHZ3orCHM4XcHPF0ccTD2QFPFwc8XRxv3zsYlzs5ZDG5q3gsSUIkhBCiQEhO03PiehzHr8ZS+NBsOoZ/DcC36R2Zdq4/nDMAV0xe4+ZkT8lCbpT0daOUrxslCrni5+VMYXdninhoiY+Pq6PUpNgASYiEEELkO2l6A0euxHL0SgzHrsVx7GosZyMS0BsMjHZYznMOvwEwK70H8x370qiUF6VuJz0lfe8kQIXdnYz9YIRtk4RICCFEvnDzVipbz0Sw6WQE285EEp+cfk8JxbuuS3hO/QnA6Wpj6d1uHK94u0jSIx5IEiIhhBCPJaUUZ8IT2HQqnM0nIzh46SZ3jyIv5OZI7VKFqFbcm2qBHjQ9Mx33I1oyRIcZVGw43DqBi3xJEiIhhBCPjZR0PXvOR7P5ZDibTkVw5WaSyfpKAZ60rexHm0r+1Crpg72dDvTpsPJ/cGQxoIOun0OdAdY5AJFvSUIkhBDC6sJik/lx9wUW77tETOKdoe5ODnY0LVuYNpX9aVPJj+I+rqYv1KfBb8O0Gad19tDja6jx1COOXhQEkhAJIYSwmiNXYvhhRyirjlwn/XZ7mL+XM20q+dO2kh9NyhXGzSmbn6q0ZFj2PJxZA3aO8NQ8qNzl0QUvChRJiIQQQjxSeoNiw4kwvt8Ryj8XbhqXNyjty+BmwbSv4q81hd1PaiIs6Qfnt4CDCzzzk3YdMiFySRIiIYQQj0R8chpL/7nM/F0XjH2DHOx0dKlZjMFNg6lewjtnG9KnazVD57eAozv0WwrBzfMucGETJCESQgiRp24kpDB7yzl+2X+ZhBRtqHwhN0f6NwziucZB+Hu55HxjSmmzT59dp9UMPfcblGqUR5ELWyIJkRBCiDyhlOK3g1d5d/UJY0fp8n4eDG4WTPdaxXF1sjd/o7u+gP3fAzro9Z0kQ8JiJCESQghhcZejE3ljxVH+PnsDgMqBXrzeoRItyhfJ/SSJx1fAhkna45D3pQO1sChJiIQQQliM3qBYsOsCH68/TWKqHicHO0a1K8+w5mVwtH+IC51e2gu/3Z5oscFwaPSSZQIW4jZJiIQQQljE6bB4Jiw/wqHLMQA0CPblg57VKVPU4+E2HHUOlvQFfQpU7AhPTge5FIewMEmIhBBCPJSUdD1fbTnHnK3/kaZXeDo78HrHSvStX+rhrxJ/KwoW9YbEKChWW+s3ZJeLvkdCPIAkREIIIXLtwMWbvL78CGcjEgBoV9mf97pXI8DbjJFj2UlL1uYaij4P3qWg71Jwcn/47QqRBUmIhBBCmE1vUHy49hTf/n0epaCIhxNTu1alU/VAy1xZ3mCA31+Ey3vA2Rv6LwNP/4ffrhDZkIRICCGEWVLTDYz+5RCrj1wHoFedErzVqTKF3J0st5NNb2ujyuwcoc9P4FfJctsWIguSEAkhhMixxNR0hi88wN9nb+Bor+PTp2vRpWYxy+5k/w+wc6b2uNuXENzCstsXIguSEAkhhMiRmMRUBs3/h38vxeDqaM/c5+rSskJRy+7kzHpYPVZ73OoNqNnHstsXIhuSEAkhhHig8LhkBny/j9Ph8Xi7OjJvUH3qlCpk2Z1En4dfB4MyQK3+0HK8ZbcvxH1IQiSEEOK+Lty4xbPf7+XKzST8PJ1ZOKQhFQM8LbuT9FQtGUqNh1KNofNMmWtIPFIPMW3ow9Pr9UyaNIng4GBcXV0pW7Ys7777LkopYxmlFJMnTyYwMBBXV1fatWvH2bNnTbYTHR1N//798fLywsfHhyFDhpCQkGBS5siRIzRv3hwXFxdKlizJjBkzHskxCiFEfnbiWhy95+7mys0kggq7sfylJpZPhgA2vwPX/gUXH22uIQcLdtAWIgesmhB9+OGHzJkzhy+//JKTJ0/y4YcfMmPGDL744gtjmRkzZvD5558zd+5c9u7di7u7OyEhISQnJxvL9O/fn+PHj7NhwwZWrVrF9u3beeGFF4zr4+LieOKJJwgKCuLAgQN89NFHTJ06lW+++eaRHq8QQuQn+0Kjeeab3dxISKFyoBfLXmxMSV83y+/ov43aRVsBus8G7xKW34cQD6BTd1fHPGKdO3fG39+f77//3risV69euLq68tNPP6GUolixYowdO5bXXnsNgNjYWPz9/Zk/fz59+vTh5MmTVKlShX/++Yd69eoBsHbtWjp27MiVK1coVqwYc+bM4c033yQsLAwnJ+2/jtdff53ff/+dU6dOPTDOuLg4vL29iY2NxcvLKw/eCSGEeLxsPhXOSz8dJCXdQP3ShfhuYH28XR0tv6P4cJjbFG5FQv1h0Oljy+9D2Cxzfr+tWkPUpEkTNm3axJkzZwA4fPgwO3bsoEOHDgCEhoYSFhZGu3btjK/x9vamYcOG7N69G4Ddu3fj4+NjTIYA2rVrh52dHXv37jWWadGihTEZAggJCeH06dPcvHkzU1wpKSnExcWZ3IQQwlb8/u9VXvjxACnpBtpU8uPHwQ3zJhkyGGDFcC0Z8qsKT7xr+X0IkUNW7VT9+uuvExcXR6VKlbC3t0ev1zNt2jT69+8PQFhYGAD+/qazk/r7+xvXhYWF4efnZ7LewcEBX19fkzLBwcGZtpGxrlAh05ES06dP5+2337bQUQohRP7xx6GrjFp6CIButYrx8VM1H+4q9fez+ws4vwUcXOGpeeDomjf7ESIHrFpD9Msvv7Bo0SIWL17MwYMHWbBgAR9//DELFiywZlhMnDiR2NhY4+3y5ctWjUcIIR6FU2FxTFh+BIBnG5Xis6dr5V0ydOUAbHpHe9zhQyhaMW/2I0QOWbWGaNy4cbz++uv06aNNvFW9enUuXrzI9OnTGThwIAEBAQCEh4cTGBhofF14eDi1atUCICAggIiICJPtpqenEx0dbXx9QEAA4eHhJmUynmeUuZuzszPOzs6WOUghhMgH4pPTeOmngySnGWhevghvd6328Feqz05yHCwfDIZ0qNId6gzIm/0IYQar1hAlJiZiZ2cagr29PQaDAYDg4GACAgLYtGmTcX1cXBx79+6lcePGADRu3JiYmBgOHDhgLLN582YMBgMNGzY0ltm+fTtpaWnGMhs2bKBixYqZmsuEEMLWKKUYt+wIoTduUczbhVl9amOfV8mQUrBqNNy8oF3BvsssmW9IPBasmhB16dKFadOmsXr1ai5cuMCKFSv49NNP6dGjBwA6nY5Ro0bx3nvvsXLlSo4ePcqAAQMoVqwY3bt3B6By5co8+eSTDBs2jH379rFz505GjhxJnz59KFZMu75Ov379cHJyYsiQIRw/fpylS5cya9YsxowZY61DF0KIx8Z3f4ey9ngYjvY6vupfB19LXqT1XocWw7FfQWcPvb8HV5+825cQ5lBWFBcXp1599VVVqlQp5eLiosqUKaPefPNNlZKSYixjMBjUpEmTlL+/v3J2dlZt27ZVp0+fNtlOVFSU6tu3r/Lw8FBeXl5q0KBBKj4+3qTM4cOHVbNmzZSzs7MqXry4+uCDD3IcZ2xsrAJUbGzswx2wEEI8Zvacu6HKTFytgiasUj/uCs3bnUWeUeq9AKWmeCm17aO83ZcQyrzfb6vOQ5RfyDxEQoiCKCIumU5f7CAyPoXutYrx2TO10OVV81V6CnzXFsKOalevf+53sLPPm30JcVu+mYdICCGEdaTpDYxYfJDI+BQq+nvyfs/qeZcMAWyYoiVDboWhxzeSDInHjiREQghhg2asPcU/F27i4ezAnGfr4OaUh4OOT/0Fe+doj7vPAa/A+5cXwgokIRJCCBvz19HrfPt3KAAfP1WDMkU98m5nESfht2Ha40YvQ4WQvNuXEA9BEiIhhLAh5yITGP+rNvniCy3K8GS1PKytSYyGn/tCagKUbg7t38m7fQnxkCQhEkIIG5GYms5LPx0gISWdBsG+jA/Jw9mh9enw6yC4GQo+peCpBWCfB9dDE8JCJCESQggboJRi4m9HOROeQFFPZ77sVxuHvLosB8D6t+D8VnB0h75LwL1w3u1LCAuQhEgIIWzAT3su8seha9jb6fiqXx38PF3ybmcHF97pRN3za/Cvmnf7EsJCJCESQogC7nJ0Iu+tPgnAxA6VaBDsm3c7u7RHuzQHQKs3oHKXvNuXEBYkCZEQQhRwH6w5RUq6gcZlCjOkWXDe7Sj2Cix9FgxpUKUbtBiXd/sSwsIkIRJCiAJs7/koVh+9jp0OJnepkneTL6YmwpJ+cCsS/Ktr8w3ZyU+MyD/k2yqEEAWU3qB4+88TAPRtUIrKgXl06SGlYOVIuH5Ym4m672Jwcs+bfQmRRyQhEkKIAmrZ/sucuB6Hp4sDY9pXyLsd7fgUji0HOwd4eqE2zF6IfEYSIiGEKIDiktP4eP1pAEa1q0BhD+e82dHpNbDpXe1xx4+gdNO82Y8QeUwSIiGEKIC+3PwfNxJSKVPUnQGNg/JmJxGnYPkwQEH9oVBvcN7sR4hHQBIiIYQoYEJv3GLeTu1aZZM6VcExLyZgTI6DJX0hNV67LMeTH1h+H0I8QpIQCSFEATNt9QnS9IpWFYvSupKf5XegFPz5CkSfB++SclkOUSBIQiSEEAXI32cj2XgyAgc7HW91qpI3O9n/AxxfoXWifmq+XJZDFAiSEAkhRAGRrjfwzu1h9s81DqKcn4fld3L9CKydqD1uNxVK1LP8PoSwAkmIhBCigFi09xJnIxIo5ObIqLZ5MMw+JR6WPQ/6FKjwJDQeafl9CGElkhAJIUQBEJOYymcbzwAw5omKeLtZuE+PUto1yqLPgVcJbSbqvJr1WggrkIRICCEKgJkbzxKTmEalAE/61i9p+R0c/BGOLgOdPfT+Adzy8AKxQliBJERCCJHPnQmPZ+GeiwBM7lwFB0sPsw8/DmvGa4/bToJSDS27fSEeA5IQCSFEPqaU4t1VJ9AbFE9U8adJuSKW3UFKgtZvKD0ZyrWHJq9advtCPCYkIRJCiHxs08kI/j57Ayd7O97sVNnyO/jrNbhxBjwDocdcuYK9KLDkmy2EEPlUarqBaX+dBGBws2CCClv4CvOHFsPhn0FnB72+B3cL1z4J8RiRhEgIIfKpn/ddIvTGLYp6OjOyTTnLbjziFKweqz1u/YZctFUUeJIQCSFEPqQ3KL7bcR6AV9qUw8PZwXIbT03U+g2lJUKZ1tBsrOW2LcRjShIiIYTIh9YdD+NydBKF3BzpXdfCw+zXjIfIk+DhDz2/kX5DwibIt1wIIfIZpRTfbNdqh55rFISrk73lNn54Cfy7ENBBz2/BIw8uDivEY0gSIiGEyGcOXLzJocsxODnY8Vzj0pbbcNhR+HOU9rjV61CmpeW2LcRjThIiIYTIZ779W6sd6lGrOEU9nS2z0aSbsPRZSE/S5htqMd4y2xUin5CESAgh8pELN26x/kQ4AEObB1tmowYDrHgRbl4An1LSb0jYJPnGCyFEPvLDzlCUgtYVi1Le39MyG/37EzizFhxc4Jmf5DplwiZJQiSEEPnEzVupLNt/BYBhzctYZqP/bYQt07THnT6FwJqW2a4Q+YwkREIIkU8s2nuRpDQ9VQK9aFy28MNv8OZFWD4UUFB3ENTu//DbFCKfMjshOnjwIEePHjU+/+OPP+jevTtvvPEGqampFg1OCCGEJiVdz4Ld2hXth7UIRqfTPdwG05Lhl+e0ztTF6kCHDy0QpRD5l9kJ0fDhwzlz5gwA58+fp0+fPri5ubFs2TLGj5dRCUIIkRf+OHSNyPgUArxc6Fyj2MNv8K/X4PphcCsMT/8IDhYarSZEPmV2QnTmzBlq1aoFwLJly2jRogWLFy9m/vz5LF++3NLxCSGEzVNK8d3tofaDmpbG0f4hezscWKBNvqizg94/gI+FZ7oWIh8y+69KKYXBYABg48aNdOzYEYCSJUty48YNy0YnhBCCbWciOROegLuTPX0alHq4jV09oNUOAbSZBGVaPXR8QhQEZidE9erV47333mPhwoVs27aNTp06ARAaGoq/v7/FAxRCCFv33d+hAPRpUApvV8fcb+hWFPwyEPSpUKkzNBttoQiFyP/MTohmzpzJwYMHGTlyJG+++SblypUD4Ndff6VJkyYWD1AIIWzZiWtx7PjvBvZ2OgY1LZ37DRn0sHwIxF4G37LQfTY8bMdsIQoQB3NfUKNGDZNRZhk++ugj7O0teIFBIYQQfLdD6zvUoVoAJQq55X5DW6fD+S3g6KZNvujibaEIhSgYzE6IMqSmphIREWHsT5ShVKmHbN8WQggBQFhsMisPXQMeciLGG//B359qj7t+Af5VLBCdEAWL2QnRmTNnGDJkCLt27TJZrpRCp9Oh1+stFpwQQtiy+bsukG5QNAj2pWZJn9xvaOv7oPRQ4Umo3tti8QlRkJidEA0aNAgHBwdWrVpFYGDgw08OJoQQIpNbKeks3nt7IsaHqR0KOwbHbk+J0vpNC0QmRMFkdkJ06NAhDhw4QKVKlfIiHiGEEMAv+y8Tl5xOmSLutK3kl/sNbXlfu6/aAwJrWCY4IQogs0eZValSReYbEkKIPJSuN/D9Dm2o/eBmwdjZ5bIm/uoBOL1am4Cx1RsWjFCIgsfshOjDDz9k/PjxbN26laioKOLi4kxuQgghHs664+FcuZmEr7sTveqUyP2GNr+n3dfoA0UrWCY4IQoos5vM2rVrB0Dbtm1NlkunaiGEsIyMofbPNgrC1SmX05lc2AnnNoOdA7SaYMHohCiYzE6ItmzZkhdxCCGEAE6HxfPvpRgc7HQ81ygodxtR6k7tUJ0BUKi0xeIToqAyOyFq2bJlXsQhhBACWH7wCgBtK/tR1DOXV6A/twku7QIHF2gxzoLRCVFw5WpixpiYGL7//ntOnjwJQNWqVRk8eDDe3jLzqRBC5Fa63sBvB68C0LtuLq9Af3ftUP2h4FXMQtEJUbCZ3al6//79lC1bls8++4zo6Giio6P59NNPKVu2LAcPHsyLGIUQwiZsPxvJjYQUCrs70api0dxt5NRquPYvOLpD01EWjU+IgszsGqLRo0fTtWtXvv32WxwctJenp6czdOhQRo0axfbt2y0epBBC2IJfD2jNZd1rF8fR3uz/V7ULuG6Zpj1u9BJ45DKpEsIGmZ0Q7d+/3yQZAnBwcGD8+PHUq1fPosEJIYStuHkrlY0nIgDoXTeXQ+2P/QYRJ7QLtzb5nwWjE6LgM/tfEC8vLy5dupRp+eXLl/H09LRIUEIIYWv+PHKNVL2BqsW8qBzoZf4G9OnaNctAS4ZcfSwanxAFndkJ0TPPPMOQIUNYunQply9f5vLlyyxZsoShQ4fSt2/fvIhRCCEKvIzmslzXDh1eDNHnwa0wNHzJgpEJYRvMbjL7+OOP0el0DBgwgPT0dAAcHR156aWX+OCDDyweoBBCFHSnw+I5ciUWR3sd3WoVN38D6SmwbYb2uNkYcPawbIBC2ACzEyInJydmzZrF9OnTOXfuHABly5bFzc3N4sEJIYQtyJh7qE0lP3zdnczfwIEFEHsZPAOh/hALRyeEbcjVPEQAbm5uVK9e3ZKxCCGEzUl72LmHUhPh74+1xy3GgaOrBaMTwnbkqA9Rz549jRdu7dmz531v5rp69SrPPvsshQsXxtXVlerVq7N//37jeqUUkydPJjAwEFdXV9q1a8fZs2dNthEdHU3//v3x8vLCx8eHIUOGkJCQYFLmyJEjNG/eHBcXF0qWLMmMGTPMjlUIISxt+5mHnHto3zeQEA4+paD2c5YPUAgbkaOEyNvbG51OB2ijzLy9vbO9mePmzZs0bdoUR0dH1qxZw4kTJ/jkk08oVKiQscyMGTP4/PPPmTt3Lnv37sXd3Z2QkBCSk5ONZfr378/x48fZsGEDq1atYvv27bzwwgvG9XFxcTzxxBMEBQVx4MABPvroI6ZOnco333xjVrxCCGFpDzX3UHIc7JypPW41ERxy0dwmhNAoK5owYYJq1qxZtusNBoMKCAhQH330kXFZTEyMcnZ2Vj///LNSSqkTJ04oQP3zzz/GMmvWrFE6nU5dvXpVKaXU7NmzVaFChVRKSorJvitWrJijOGNjYxWgYmNjzTo+IYS4n+iEFFXujdUqaMIqdeJaLs4va15XaoqXUl/UU0qfbvkAhcjnzPn9NnvYfZs2bYiJicm0PC4ujjZt2pi1rZUrV1KvXj2eeuop/Pz8qF27Nt9++61xfWhoKGFhYbRr1864zNvbm4YNG7J7924Adu/ejY+Pj8mkkO3atcPOzo69e/cay7Ro0QInpzv/PYWEhHD69Glu3ryZKa6UlBTi4uJMbkIIYWkrD18jTa9yN/fQ+W2wZ7b2OOR9sLO3fIBC2BCzE6KtW7eSmpqaaXlycjJ///23Wds6f/48c+bMoXz58qxbt46XXnqJV155hQULFgAQFhYGgL+/v8nr/P39jevCwsLw8/MzWe/g4ICvr69Jmay2cfc+7jZ9+nSTZsCSJXN5kUUhhLiPXM89lBwLv7+sPa77PJRvb9nAhLBBOR5lduTIEePjEydOmCQSer2etWvXUry4efNnGAwG6tWrx/vva7Or1q5dm2PHjjF37lwGDhxo1rYsaeLEiYwZM8b4PC4uTpIiIYRFnQqL4+jVXM49tGYCxF2BQsHwxLS8CVAIG5PjhKhWrVrodDp0Ol2WTWOurq588cUXZu08MDCQKlWqmCyrXLkyy5cvByAgIACA8PBwAgMDjWXCw8OpVauWsUxERITJNtLT04mOjja+PiAggPDwcJMyGc8zytzN2dkZZ2dns45FCCHMsfxALuceOrESDv8MOjvoMVcmYRTCQnLcZBYaGsq5c+dQSrFv3z5CQ0ONt6tXrxIXF8fgwYPN2nnTpk05ffq0ybIzZ84QFBQEQHBwMAEBAWzatMm4Pi4ujr1799K4cWMAGjduTExMDAcOHDCW2bx5MwaDgYYNGxrLbN++nbS0NGOZDRs2ULFiRZMRbUII8Sik6Q2s+PcaYObcQ/HhsGqU9rjpq1CqkeWDE8JW5X0f7+zt27dPOTg4qGnTpqmzZ8+qRYsWKTc3N/XTTz8Zy3zwwQfKx8dH/fHHH+rIkSOqW7duKjg4WCUlJRnLPPnkk6p27dpq7969aseOHap8+fKqb9++xvUxMTHK399fPffcc+rYsWNqyZIlys3NTX399dc5ilNGmQkhLGnjiTAVNGGVqvPOepWars/ZiwwGpRY9rY0qm91UqbSUB79GCBtnzu+32QnR+++/r77//vtMy7///nv1wQcfmLs59eeff6pq1aopZ2dnValSJfXNN9+YrDcYDGrSpEnK399fOTs7q7Zt26rTp0+blImKilJ9+/ZVHh4eysvLSw0aNEjFx8eblDl8+LBq1qyZcnZ2VsWLFzcrVkmIhBCW9OLC/Spowir1zp/Hc/6iAwu0ZOidIkqFHcu74IQoQMz5/dYppZQ5NUqlS5dm8eLFNGnSxGT53r176dOnD6GhoRarvXpcxMXF4e3tTWxsLF5eZg6NFUKIu9y8lUqD9zeSplesebV5zobbR4fC3GaQmgDt39Gay4QQD2TO77fZw+7DwsJMOjhnKFq0KNevXzd3c0IIYVPMnnvIoIffX9KSoVJNoPHIvA9SCBtkdkJUsmRJdu7cmWn5zp07KVasmEWCEkKIgsrsuYd2fwmXdoOTB/SYIxMwCpFHzL7a/bBhwxg1ahRpaWnG4febNm1i/PjxjB071uIBCiFEQWH23ENhx2Dze9rjJ6dDodJ5Gp8QtszshGjcuHFERUXx8ssvG2esdnFxYcKECUycONHiAQohREFh1txD6SmwYjjoU6FCB7mSvRB5zOyESKfT8eGHHzJp0iROnjyJq6sr5cuXl4kMhRDiPsyee2jrdAg/Bm6FoevnoNPlcYRC2Daz+xBlCAsLIzo6mrJly+Ls7IyZg9WEEMKmbD8TyY2EFAq7O9GqYtH7F760B3bO0h53mQUefvcvL4R4aGYnRFFRUbRt25YKFSrQsWNH48iyIUOGSB8iIYTIxp+HtdqhLjWL4Wh/n1Nv6i1Y8SIoA9TsB5W7PKIIhbBtZidEo0ePxtHRkUuXLuHm5mZc/swzz7B27VqLBieEEAVBcpqeDSe06yd2qfmA0bhb3oeboeBVHDp88AiiE0JALvoQrV+/nnXr1lGihOmQ0fLly3Px4kWLBSaEEAXFtjOR3ErVU8zbhdolfbIveGU/7JmtPe4yC1y8H0l8Qohc1BDdunXLpGYoQ3R0tHSsFkKILKw6onUt6Fg9EDu7bDpHp6fAHyO0prIafaB8+0cYoRDC7ISoefPm/Pjjj8bnOp0Og8HAjBkzaN26tUWDE0KI/C4pVc+mk1pzWacamWf5N9r+MUSeAvei2pxDQohHyuwmsxkzZtC2bVv2799Pamoq48eP5/jx40RHR2c5g7UQQtiyracjSEzVU9zHlVrZNZeFHYUdn2qPO34Mbr6PLD4hhMbsGqJq1apx5swZmjVrRrdu3bh16xY9e/bk33//pWzZsnkRoxBC5FurjmrNZZ1rBKLLai4hfTr8MRIM6VCpM1Tp9ogjFEJADmuIevbsyfz58/Hy8uLHH3/kmWee4c0338zr2IQQIl9LTE1n88kI4D7NZbu/gOuHtA7UnT6RCRiFsJIc1RCtWrWKW7duATBo0CBiY2PzNCghhCgItpyKJClNT0lfV6oXz2LE2I2zsOV2f6GQ6eAZ8GgDFEIY5aiGqFKlSkycOJHWrVujlOKXX37By8sry7IDBgywaIBCCJFfrTqiTcbYqXqxzM1lBgOs/B/oU6BsW6jVzwoRCiEy5Cghmjt3LmPGjGH16tXodDreeuutLNvCdTqdJERCCAHcSkln8ymtuaxzVs1l/3wHl3aDkwd0mSlNZUJYWY4SoiZNmrBnzx4A7OzsOHPmDH5+cm0dIYTIzqZTEaSkGyhd2I2qxe6pUY+5BBunao/bTQWfUo86PCHEPcweZRYaGkrRog+4MKEQQti41RnNZfeOLlMK/nwV0m5BqSZQb4iVIhRC3M3seYiCgoKIiYlh3759REREYDAYTNZLk5kQwtYlpKSz5XQkoPUfMnFoMZzbDA4u0PULsDP7/1IhRB4wOyH6888/6d+/PwkJCXh5eZn85yN9iIQQAjaeCCc13UCZIu5UDvS8syI+DNZN1B63mghFylknQCFEJmb/azJ27FgGDx5MQkICMTEx3Lx503iLjo7OixiFECJfybh2mclkjErB6rGQHAuBtaDxSOsFKITIxOyE6OrVq7zyyitZXuBVCCFsXVxyGtvP3G4uq3FXc9nJlXBqFdg5QLevwN7sCnohRB4yOyEKCQlh//79eRGLEELkextPhJOqN1DOz4MK/h7awuQ4WDNBe9xsNARUs16AQogsmf0vSqdOnRg3bhwnTpygevXqODo6mqzv2rWrxYITQoj8ZvXt5rJO1e9qLtvyPsRfh0LB0Pw1K0YnhMiO2QnRsGHDAHjnnXcyrdPpdOj1+oePSggh8qHYxDS2n9Way4yTMV47BPu+1h53+gQcXawTnBDivsxOiO4dZi+EEEKz/kQYaXpFRX9Pyvt7gkEPq0aDMkC1XlCurbVDFEJkQybAEEIIC1l99HZzWUbt0P4f4NpBcPbWLt4qhHhs5biG6PPPP89RuVdeeSXXwQghRH4Vk5jKjrM3AOhYPVCbc2jT7a4FbSeBp78VoxNCPEiOE6LPPvvsgWV0Op0kREIIm7T+eDjpBkWlAE/K+XnAr69AShwUqwP1Bls7PCHEA+Q4IQoNDc3LOIQQIl9bdbu5rEvNYvDfJji2HHR20PkzsLO3cnRCiAeRmcGEEOIhRd9KZed/WnNZp8o+sLS3tqLBcChWy2pxCSFyTjpVCyHEQ1p3PAy9QVG1mBelT34DN0PBMxBav2Ht0IQQOSQJkRBCPKSMyRj7l0uBHbf7Wz75Abh4WTEqIYQ5JCESQoiHEJWQwq5zNwBFz2ufgT4Vyj8BVbpZOzQhhBkkIRJCiIew9ngYBgWvFDmIy5Ud4OACHT+CjMt2CCHyhVwlROfOneOtt96ib9++REREALBmzRqOHz9u0eCEEOJxt+rwdbxI4KWUH7QFLcdDodJWjUkIYT6zE6Jt27ZRvXp19u7dy2+//UZCQgIAhw8fZsqUKRYPUAghHlc3ElLYGxrFBIeluKbdhCIVofH/rB2WECIXzE6IXn/9dd577z02bNiAk5OTcXmbNm3Ys2ePRYMTQojH2frj4dTiDP0dNmkLOn8GDk73f5EQ4rFkdkJ09OhRevTokWm5n58fN27csEhQQgiRH6w7eoX3HOdpT2r1h9JNrRuQECLXzE6IfHx8uH79eqbl//77L8WLF7dIUEII8bi7eSuV4AtLqWJ3Eb2zD7R/x9ohCSEegtkJUZ8+fZgwYQJhYWHodDoMBgM7d+7ktddeY8CAAXkRoxBCPHa2/3uC0fa/AGDfbhK4F7FyREKIh2F2QvT+++9TqVIlSpYsSUJCAlWqVKFFixY0adKEt956Ky9iFEKIx47Pnul46xKJ8KgIdQdZOxwhxEPSKaVUbl546dIljh07RkJCArVr16Z8+fKWju2xERcXh7e3N7GxsXh5ycyzQti6hHN78Vj4BABXev5OiRqtrRyRECIr5vx+m31x1x07dtCsWTNKlSpFqVKlch2kEELkSwYDqStHA7DesQ1PSDIkRIFgdpNZmzZtCA4O5o033uDEiRN5EZMQQjy+/v0R39jjxClXLtQeZ+1ohBAWYnZCdO3aNcaOHcu2bduoVq0atWrV4qOPPuLKlSt5EZ8QQjw+EqMxbHwbgJnpvWlVt7qVAxJCWIrZCVGRIkUYOXIkO3fu5Ny5czz11FMsWLCA0qVL06ZNm7yIUQghHg9bpmGXFM1pQwl2FOpOeT8Pa0ckhLCQh7q4a3BwMK+//joffPAB1atXZ9u2bZaKSwghHi/XD8N+7XplU9MHElKjJDq5gKsQBUauE6KdO3fy8ssvExgYSL9+/ahWrRqrV6+2ZGxCCPF4UAr+GgfKwGpDY3YbqtKhWqC1oxJCWJDZo8wmTpzIkiVLuHbtGu3bt2fWrFl069YNNze3vIhPCCGs78hSuLyXdHs33k3uR+nCblQO9LR2VEIICzI7Idq+fTvjxo3j6aefpkgRmZlVCFHAJcfB+kkArPJ5lrBbhXmpeqA0lwlRwJidEO3cuTMv4hBCiMfTtg/hVgQG33JMjmgBQEdpLhOiwMlRQrRy5Uo6dOiAo6MjK1euvG/Zrl27WiQwIYSwuoiTsGcOAAerTCBuox0lCrlSrbjMWC9EQZOjhKh79+6EhYXh5+dH9+7dsy2n0+nQ6/WWik0IIazH2JFaD5U682NkeeAaHaoFSHOZEAVQjhIig8GQ5WMhhCiwjq+AC3+Dgwspbd9l0xdnAOhQXZrLhCiIzB52/+OPP5KSkpJpeWpqKj/++KNFghJCCKtKSYD1b2mPm43m70h3bqXqCfR2oVYJH6uGJoTIG2YnRIMGDSI2NjbT8vj4eAYNGmSRoIQQwqq2Toe4q+BTCpq+yl/HrgPwZLUA7OykuUyIgsjshEgplWX7+ZUrV/D29rZIUEIIYTVX9sOe2drjjp+QqnNmw4lw7ak0lwlRYOV42H3t2rXR6XTodDratm2Lg8Odl+r1ekJDQ3nyySfzJEghhHgk0lPgjxGgDFDjGajwBDtPRxCfnI6fpzN1SxWydoRCiDyS4xqi7t27061bN5RShISE0K1bN+OtT58+fP311/z000+5DuSDDz5Ap9MxatQo47Lk5GRGjBhB4cKF8fDwoFevXoSHh5u87tKlS3Tq1Ak3Nzf8/PwYN24c6enpJmW2bt1KnTp1cHZ2ply5csyfPz/XcQohCrC/P4HIU+BWBJ78AIA1R6W5TAhbkOMaoilTpgBQunRpnnnmGVxcXCwWxD///MPXX39NjRo1TJaPHj2a1atXs2zZMry9vRk5ciQ9e/Y0Tg6p1+vp1KkTAQEB7Nq1i+vXrzNgwAAcHR15//33AQgNDaVTp068+OKLLFq0iE2bNjF06FACAwMJCQmx2DEIIfK5sGNaQgTQ8SNw8yVNb2D97eYyuXaZEAWcsrL4+HhVvnx5tWHDBtWyZUv16quvKqWUiomJUY6OjmrZsmXGsidPnlSA2r17t1JKqb/++kvZ2dmpsLAwY5k5c+YoLy8vlZKSopRSavz48apq1aom+3zmmWdUSEhIjmOMjY1VgIqNjc3tYQohHmfpaUp93VKpKV5K/dxPKYNBKaXU9jMRKmjCKlXnnfUqXW+wboxCCLOZ8/ttdqdqvV7Pxx9/TIMGDQgICMDX19fkZq4RI0bQqVMn2rVrZ7L8wIEDpKWlmSyvVKkSpUqVYvfu3QDs3r2b6tWr4+/vbywTEhJCXFwcx48fN5a5d9shISHGbWQlJSWFuLg4k5sQogDbMxuu/QvO3tDpE7g9cOSvo2EAhFQLwF6ay4Qo0MxOiN5++20+/fRTnnnmGWJjYxkzZgw9e/bEzs6OqVOnmrWtJUuWcPDgQaZPn55pXVhYGE5OTvj4+Jgs9/f3JywszFjm7mQoY33GuvuViYuLIykpKcu4pk+fjre3t/FWsmRJs45LCJGPRJ2DLdO0xyHTwDMAgHS9gfXHtfOIXLtMiILP7IRo0aJFfPvtt4wdOxYHBwf69u3Ld999x+TJk9mzZ0+Ot3P58mVeffVVFi1aZNH+SJYwceJEYmNjjbfLly9bOyQhRF4wGGDlK5CeDGVaQe1njav2XYgm6lYqhdwcaVjG/NpvIUT+YnZCFBYWRvXq1QHw8PAwTtLYuXNnVq9enePtHDhwgIiICOrUqYODgwMODg5s27aNzz//HAcHB/z9/UlNTSUmJsbkdeHh4QQEaP/BBQQEZBp1lvH8QWW8vLxwdXXNMjZnZ2e8vLxMbkKIAujAPLi4AxzdoMssY1MZwJrbzWVPVAnA0d7sU6UQIp8x+6+8RIkSXL+uDUMtW7Ys69evB7SRYs7OzjneTtu2bTl69CiHDh0y3urVq0f//v2Njx0dHdm0aZPxNadPn+bSpUs0btwYgMaNG3P06FEiIiKMZTZs2ICXlxdVqlQxlrl7GxllMrYhhLBRsVdggzZ6lrZToFBp4yq9QbH2dnNZh+oBVghOCPGo5XjYfYYePXqwadMmGjZsyP/+9z+effZZvv/+ey5dusTo0aNzvB1PT0+qVatmsszd3Z3ChQsblw8ZMoQxY8bg6+uLl5cX//vf/2jcuDGNGjUC4IknnqBKlSo899xzzJgxg7CwMN566y1GjBhhTM5efPFFvvzyS8aPH8/gwYPZvHkzv/zyi1m1WUKIAkYpWDUaUuOhRANoMMxk9YGLN4mMT8HLxYEmZYtYKUghxKNkdkL0wQcfGB8/88wzxlFf5cuXp0uXLhYN7rPPPsPOzo5evXqRkpJCSEgIs2fPNq63t7dn1apVvPTSSzRu3Bh3d3cGDhzIO++8YywTHBzM6tWrGT16NLNmzaJEiRJ89913MgeRELbs6DI4ux7snaDbl2Bnb7L6z8PXAGhfJQAnB2kuE8IW6JRSytpBPO7i4uLw9vYmNjZW+hMJkd8lRMJXDSApGtpMghavmaxO1xto+P4mom6lMn9QfVpV9LNSoEKIh2XO73eOaohWrlyZ45137do1x2WFEOKRWzNOS4YCqkPTVzOt3nUuiqhbqfi6O9G0nDSXCWErcpQQde/ePUcb0+l06PX6h4lHCCHyzslVcHwF6Oyh65dg75ipyMrbzWUdq8voMiFsSY4SIoPBkNdxCCFE3kqKgdVjtcdNX4FitTIVSU7Ts+6YNrqsa83ijy42IYTVyb8/QgjbsGESJIRB4XLQ8vUsi2w9HUl8SjqB3i7UCyr0iAMUQliT2aPM7h7BlZXJkyfnOhghhMgT57fCwR+1x12/BMesZ8f/84jWXNa5RiB2cu0yIWyK2QnRihUrTJ6npaURGhqKg4MDZcuWlYRICPF4Sb2lXZ4DoP4wCMp6UtZbKelsOqnNai/NZULYHrMTon///TfTsri4OJ5//nl69OhhkaCEEMJiNk+DmIvgXRLaTcm22IYT4SSnGQgu4k614jK9hhC2xiJ9iLy8vHj77beZNGmSJTYnhBCWcfkf2HN7MtfOM8HZM9uiGaPLutQshk4nzWVC2BqLdarOuDK8EEI8FtJTYOVIQEHNvlC+XbZFb95KZfuZSAC61gx8RAEKIR4nZjeZff755ybPlVJcv36dhQsX0qFDB4sFJoQQD+XvTyHyFLgXhZD371t07fEw0g2KyoFelPPLvhZJCFFwmZ0QffbZZybP7ezsKFq0KAMHDmTixIkWC0wIIXIt/Dj8/Yn2uONH4OZ73+IrD2nNZV1rFsvryIQQjymzE6LQ0NC8iEMIISxDnw5/jARDGlTqDFW637d4eFwye0KjAOgizWVC2CyZmFEIUbDsnQPXDoKzN3T8GB7QQXrVkesoBXWDClGikNsjClII8bgxu4YoOTmZL774gi1bthAREZHpsh4HDx60WHBCCGGWqHPaMHuAkPfA68E1PsbRZTWkdkgIW2Z2QjRkyBDWr19P7969adCggQxPFUI8HpSCP1+F9CQIbgm1n3vgSy5FJXL4cgx2OuhUQ/oPCWHLzE6IVq1axV9//UXTpk3zIh4hhMidgwvgwt/g6AZdZj2wqQzuXKqjSdkiFPV0zusIhRCPMbP7EBUvXhxPTxmWKoR4jMRdg/W3J4ZtMwl8g3P0MhldJoTIYHZC9MknnzBhwgQuXryYF/EIIYR5lIJVYyAlDorXg4bDc/Sy02HxnA6Px9FeR0i1gDwOUgjxuDO7yaxevXokJydTpkwZ3NzccHR0NFkfHR1tseCEEOKBDi6AM2vAzhG6fQl29jl62crDVwFoWcEPb1fHB5QWQhR0ZidEffv25erVq7z//vv4+/tLp2ohhPVc3AWrX9Met34D/Crn6GVKKf48fB2ArrWkuUwIkYuEaNeuXezevZuaNWvmRTxCCJEzMZdg6XPaBIxVe0Cz0Tl+6eErsVyKTsTV0Z52lf3yMEghRH5hdh+iSpUqkZSUlBexCCFEzqQkwM99IfEGBNSAbrNzNKosQ0Zn6vZV/HFzMvv/QiFEAWR2QvTBBx8wduxYtm7dSlRUFHFxcSY3IYTIUwYD/P4ihB8Ddz/o+zM45XyGab1BseqIjC4TQpgy+1+jJ598EoC2bduaLFdKodPp0Ov1lolMCCGysu0DOPkn2DtBn0XgXcKsl+8NjSIiPgUvFweaVyiSR0EKIfIbsxOiLVu25EUcQgjxYMdXwLYPtcddZkHJBmZvIqMzdYdqgTg75GxEmhCi4DM7IWrZsmVexCGEEPd3/TCseEl73Hgk1Opn9iZS0w2sOSajy4QQmZmdEG3fvv2+61u0aJHrYIQQIksJEfBzP+06ZeXaQft3crWZHf9FEpOYRhEPZxqVKWzhIIUQ+ZnZCVGrVq0yLbt7LiLpQySEsKj0FFj6LMRdgcLlodf3OZ588V4Zo8s61wjE3k7mUBNC3GH2KLObN2+a3CIiIli7di3169dn/fr1eRGjEMJWZVyW4/JecPaGvkvA1SdXm4pPTmPd8XAAusjoMiHEPcyuIfL29s60rH379jg5OTFmzBgOHDhgkcCEEII9c+DQT6Czg6fmQZFyud7UH4eukZSmp5yfB3VK+VguRiFEgWB2DVF2/P39OX36tKU2J4Swdf9tgvVvao+fmAbl2t6//AMs+ecSAH3ql5RLDgkhMjG7hujIkSMmz5VSXL9+nQ8++IBatWpZKi4hhC278R8sGwTKALWfhUYvPdTmjl6J5djVOJzs7ehZx7x5i4QQtsHshKhWrVrodDqUUibLGzVqxA8//GCxwIQQNio5Fn7uAymxULIhdPrUrMtyZOXn27VDT1YLwNfdyRJRCiEKGLMTotDQUJPndnZ2FC1aFBcXF4sFJYSwUQY9/DoEos6CV3F45idwcH6oTd5KSeePf68C0KdBSUtEKYQogMxOiIKCgvIiDiGEgE1vw38bwMEV+iwGj4e/Ev2qI9e4laqndGE3GsvcQ0KIbOS4U/XmzZupUqVKlhdwjY2NpWrVqvz9998WDU4IYUOO/AI7Z2mPu38FxWpZZLM/77sMQJ8GpaQztRAiWzlOiGbOnMmwYcPw8vLKtM7b25vhw4fz6aefWjQ4IYSNuHoA/hipPW42Bqr1sshmT16P49DlGBzsdPSSztRCiPvIcUJ0+PBh45Xus/LEE0/IHERCCPPFh8GS/qBPgQodoM0ki216yT6tM/UTVf0p6vlwfZGEEAVbjhOi8PBwHB0ds13v4OBAZGSkRYISQtiItGQtGYq/DkUrQc9vwM4y06Mlper5LaMzdf1SFtmmEKLgyvGZp3jx4hw7dizb9UeOHCEwMNAiQQkhbIBSsGo0XN0PLj5aJ2qXzE3yufXX0evEJ6dTopArzcoVsdh2hRAFU44Too4dOzJp0iSSk5MzrUtKSmLKlCl07tzZosEJIQqw3V/B4cWgs4en5kPhshbd/M/77sxMbScXchVCPECOh92/9dZb/Pbbb1SoUIGRI0dSsWJFAE6dOsVXX32FXq/nzTffzLNAhRAFyH8bYcPtvkIh06Bsa4tu/mx4PPsv3sTeTsdT9WTuISHEg+U4IfL392fXrl289NJLTJw40ThTtU6nIyQkhK+++gp/f/88C1QIUUDc+A+WDb5zWY6GL1p8F0v+0Ybat6nkh7+XTBorhHgwsyZmDAoK4q+//uLmzZv8999/KKUoX748hQoVyqv4hBAFSR5cliPTLtL0LD94BYC+MjO1ECKHzJ6pGqBQoULUr1/f0rEIIQoypeD3l+9cluPphQ99WY6srDseRkxiGoHeLrSs8PAzXQshbINlxrcKIcSD7JkNp1aBnaOWDHnmTRN7Rmfqp+uVxF46UwshckgSIiFE3ru0BzZM1h4/OR1K1M2T3ZyPTGDP+WjsdPB0fWkuE0LknCREQoi8lRAJywaBIV27JEf9oXm2q6W3O1O3rFCU4j6uebYfIUTBIwmRECLvGPTw21CIvwZFKkCXWRbvRJ0hNd3Arwe0ztR9GsjM1EII80hCJITIO9s+hPNbwdENnv4RnD3zbFcbT4YTdSsVP09n2lSSztRCCPNIQiSEyBtnN8K2GdrjLrPAr3Ke7i6jM/VT9UrgaC+nNiGEeeSsIYSwvNgr8NswQEG9wVDj6Tzd3eXoRP4+ewOQC7kKIXJHEiIhhGWlp8Ky5yEpGgJrQcj0PN/lkn+02qHm5YtQ0tctz/cnhCh4JCESQljWhslw5R9w8YanF4Bj3l46I01vYNn+252ppXZICJFLkhAJISzn+ArYO0d73H0uFCqd57vceCKciPgUCrs70b6KXE9RCJE7khAJISzjxln4Y6T2uOkoqNQxz3eplOKrrf8B0LdBKZwc5JQmhMgdOXsIIR5eaiL8MgBSEyCoKbSZ9Eh2u+1MJMeuxuHqaM+gpqUfyT6FEAWTJERCiIejFKweCxEnwN0Pev8A9rm6brSZu1V8sVmrHerfsBSFPSx/oVghhO2QhEgI8XAOzIPDi0FnpyVDngGPZLd7zkdz4OJNnBzsGNaizCPZpxCi4JKESAiRe1cPwJoJ2uO2UyC4+SPb9ZdbzgLwTL2S+Hvl7Ug2IUTBZ9WEaPr06dSvXx9PT0/8/Pzo3r07p0+fNimTnJzMiBEjKFy4MB4eHvTq1Yvw8HCTMpcuXaJTp064ubnh5+fHuHHjSE9PNymzdetW6tSpg7OzM+XKlWP+/Pl5fXhCFGy3omDpANCnQqXO0PTVR7brAxdvsvO/KBzsdAxvKbVDQoiHZ9WEaNu2bYwYMYI9e/awYcMG0tLSeOKJJ7h165axzOjRo/nzzz9ZtmwZ27Zt49q1a/Ts2dO4Xq/X06lTJ1JTU9m1axcLFixg/vz5TJ482VgmNDSUTp060bp1aw4dOsSoUaMYOnQo69ate6THK0SBYdDD8iEQdwV8y0L32Xl20dasfLVF6zvUs05xShSSiRiFEA9Pp5RS1g4iQ2RkJH5+fmzbto0WLVoQGxtL0aJFWbx4Mb179wbg1KlTVK5cmd27d9OoUSPWrFlD586duXbtGv7+2hwkc+fOZcKECURGRuLk5MSECRNYvXo1x44dM+6rT58+xMTEsHbt2gfGFRcXh7e3N7GxsXh5eeXNwQuRn2x+D7Z/pF20degm8K/yyHZ97Gosnb/YgZ0ONo1tRXAR90e2byFE/mLO7/dj1YcoNjYWAF9fXwAOHDhAWloa7dq1M5apVKkSpUqVYvfu3QDs3r2b6tWrG5MhgJCQEOLi4jh+/LixzN3byCiTsY17paSkEBcXZ3ITQtx2eq2WDIF20dZHmAwBzL4971DnGsUkGRJCWMxjkxAZDAZGjRpF06ZNqVatGgBhYWE4OTnh4+NjUtbf35+wsDBjmbuToYz1GevuVyYuLo6kpKRMsUyfPh1vb2/jrWTJkhY5RiHyvehQWPGC9rj+sDy/aOu9zobHs+aY9nc9onW5R7pvIUTB9tgkRCNGjODYsWMsWbLE2qEwceJEYmNjjbfLly9bOyQhrC8tCX55DpJjoUR9CHn/kYcwe+s5lIKQqv5UDPB85PsXQhRceT97Wg6MHDmSVatWsX37dkqUKGFcHhAQQGpqKjExMSa1ROHh4QQEBBjL7Nu3z2R7GaPQ7i5z78i08PBwvLy8cHV1zRSPs7Mzzs4yyZsQRhmTL4YdBbci8NQCcHB6pCFcjLrFH4euAjCydflHum8hRMFn1RoipRQjR45kxYoVbN68meDgYJP1devWxdHRkU2bNhmXnT59mkuXLtG4cWMAGjduzNGjR4mIiDCW2bBhA15eXlSpUsVY5u5tZJTJ2IYQ4gEOzIdDi+5Mvuhd/JGHMGfrOQwKWlUsSvUS3o98/0KIgs2qNUQjRoxg8eLF/PHHH3h6ehr7/Hh7e+Pq6oq3tzdDhgxhzJgx+Pr64uXlxf/+9z8aN25Mo0aNAHjiiSeoUqUKzz33HDNmzCAsLIy33nqLESNGGGt5XnzxRb788kvGjx/P4MGD2bx5M7/88gurV6+22rELkW9cPQBrxmuP20yCMi0feQjXYpJYfvAKAP9rI32HhBCWZ9Uaojlz5hAbG0urVq0IDAw03pYuXWos89lnn9G5c2d69epFixYtCAgI4LfffjOut7e3Z9WqVdjb29O4cWOeffZZBgwYwDvvvGMsExwczOrVq9mwYQM1a9bkk08+4bvvviMkJOSRHq8Q+c6tKPhloDb5YsVO0Gy0VcL4Zvt50vSKRmV8qRvka5UYhBAF22M1D9HjSuYhEjZJnw6Ln4Jzm8G3DLywFVwefVNVRHwyzT/cQkq6gUVDG9K0XJFHHoMQIn/Kt/MQCSEeEwmRsLC7lgw5uMLTC62SDAF8/3coKekGapfyoUnZwlaJQQhR8D0Wo8yEEI+RK/vhlwEQdxUc3aHXdxBQzSqh3LyVysI9FwGt75DuEV4eRAhhWyQhEkJolIL9P2hXrzekQeHy8MxP4FfJaiHN2xlKYqqeKoFetK7oZ7U4hBAFnyREQght0sVVY+DwYu155S7QbTa4WK/PXFxyGvN3XQCkdkgIkfckIRLC1kWHajNQhx3V5hlqNxWavPJIr16flYW7LxKXnE45Pw9CqgZYNRYhRMEnCZEQtuzsBlg+FJJjtBmoe/9glXmG7hURn8zX284BMKJ1WezspHZICJG3JCESwhYZDLB9Bmz9AFBQvB48vQC8SzzwpY/Ce6tOEpecTvXi3nSt+ehnxRZC2B5JiISwNUk34bcX4Ox67Xm9IfDkdHB4PK7ft/1MJCsPX8NOB+/3qI691A4JIR4BSYiEsCUXd2nJUOxlcHCBzp9BrX7WjsooOU3PW78fA2Bgk9JyzTIhxCMjCZEQtkCfBts+hL8/AWWAQqXh6R8hsKa1IzPxxeazXIpOJMDLhbFPVLR2OEIIGyIJkRAFXfR5WD4Mru7XntfsBx1ngLOndeO6x5nweL7Zfh6At7tVxcNZTk9CiEdHzjhCFFRKweEl8NdrkJoAzt7Q5TOo1svakWViMCjeXHGUNL2iXWV/GWYvhHjkJCESoiBKioFVo+H4b9rzUk2g5zfgU9KqYWXnl/2X+efCTdyc7Hm7W1VrhyOEsEGSEAlR0NzdcVpnD60nQrMxYGdv7ciydCMhhelrTgEwpn0Fivu4WjkiIYQtkoRIiIJCn3674/THdzpO9/oeStSzdmT3NW31SWKT0qgS6MXzTUpbOxwhhI2ShEiIgiDiFKwcCVf+0Z4/ph2n77Xj7A1W/HsVnQ6m96yOg72dtUMSQtgoSYiEyM/SkrUaoR0ztSvUP8Ydp++lzTl0FIABjYKoWdLHugEJIWyaJERC5Ffnt2odp6O1oepU6ACdPn5sLr/xILO3/MeFqET8vZwZGyJzDgkhrEsSIiHym1s3YN2bcGSJ9twzEDrMgMpdrH6F+pz6LyKeObcv3jq1S1W8XBytHJEQwtZJQiREfqEU/PsTbJikXY8MHTQYBm0mgYuXtaPLMaUUb6w4Rppe0aaSH09WkzmHhBDWJwmREPlB5BlYNQou7tSe+1eHLrOgRF2rhpUbyw5cYV9oNK6O9rzdtSq6fFKrJYQo2CQhEuJxlpYMOz6DHZ+CPhUc3aDVRGj0Etjnv2amyPgU3v/rJACj25enpK+blSMSQgiNJERCPI6UgtNrYN0bcDNUW1auPXT6BAoFWTe2XEpJ1/PSTweISUyjUoAng5oGWzskIYQwkoRIiMdN5GlY+zqc26w99wiAJ6dD1R75ptP0vZRSTFx+lP0Xb+Lp4sCX/ergKHMOCSEeI5IQCfG4SIqBbTNg39dgSAd7J2j0MrR47bGfYPFB5mw7x2//XsXeTsfs/nUo5+dh7ZCEEMKEJERCWJtBr40e2/QOJN7QllXoACHToHBZ68ZmAWuPhTFj7WkApnapQvPyRa0ckRBCZCYJkRDWdGkPrBkP1w9rz4tU0JrHyrWzblwWcuxqLKOXHgJgYOMgnmtc2qrxCCFEdiQhEsIa4q7BhslwdJn23NkLWr0ODV7Il6PHshIRl8zQBftJStPTokJRJnWuYu2QhECv15OWlmbtMIQFOTk5YWf38H0SJSES4lEKOwr7voEjyyA9CdBBneegzWTwKDhNSclpeob9uJ+wuGTK+XnwZb/acuFWYVVKKcLCwoiJibF2KMLC7OzsCA4OxsnJ6aG2IwmREHlNnwYnV8K+b+HS7jvLSzaCDh9AsdrWiy0PGAyKscsOc/hKLIXcHPl+YD25NIewuoxkyM/PDzc3N5kQtIAwGAxcu3aN69evU6pUqYf6XCUhEiKvxIfDgfmw/wdICNOW2TlA5a7QcDiUbJhvh9Hfz6xNZ1l95DqO9jrmPluXoMLu1g5J2Di9Xm9MhgoXLmztcISFFS1alGvXrpGeno6jY+7/+ZKESAhLUgqu/KM1ix3/HQy3+yp4+EPdQVD3efAKtGaEeWrl4WvM2nQWgGk9qtOwjPz4COvL6DPk5iYzoxdEGU1ler1eEiIhrC49FY4th71z7owYAyjRQKsNqtwVHB6ufftx9++lm7y2TDv24S3K8HS9klaOSAhT0kxWMFnqc5VejkI8jMRo+PsTmFkdfn9RS4bsnaFWf3hhKwzdANV7F/hk6FpMEsN+PEBquoF2lf0Z/2Qla4ckhM3Q6XT8/vvveb6f0qVLM3PmzDzfT1bmz5+Pj49Pnu5DEiIhciPqHKx+DT6rqk2omBCmXWKjzSQYcxK6zy5wnaWzczHqFs9+t5cbCSlUCvBkVp9a2NvJf+JCWEpkZCQvvfQSpUqVwtnZmYCAAEJCQti5cycA169fp0OHDlaOMrNHkcRYkjSZCZFTSmmjxHZ/BadWA0pb7l8dGo+Aar0KfE3QvfZfiOaFhQeIvpVKcR9Xvn++Pu7OcloRwpJ69epFamoqCxYsoEyZMoSHh7Np0yaioqIACAgIsHKEBYPUEAnxIPo0OPorfNsa5nWAU6sABeWfgAEr4cW/oVZfm0uGVh6+Rr/v9hJ9K5UaJbxZ8XITivu4WjssIQqUmJgY/v77bz788ENat25NUFAQDRo0YOLEiXTt2hUwbTK7cOECOp2OX375hebNm+Pq6kr9+vU5c+YM//zzD/Xq1cPDw4MOHToQGRlp3E+rVq0YNWqUyb67d+/O888/n21sn376KdWrV8fd3Z2SJUvy8ssvk5CQAMDWrVsZNGgQsbGx6HQ6dDodU6dOBSAlJYXXXnuN4sWL4+7uTsOGDdm6davJtufPn0+pUqVwc3OjR48exuQvL8m/ckJkJy0JDi6EXZ9D7GVtmb0z1Oyj1QgVrWjd+KxEKcVXW/7j4/VnAHiiij8z+9TCzUlOJyL/UEqRlKa3yr5dHe1z3BHYw8MDDw8Pfv/9dxo1aoSzs3OOXjdlyhRmzpxJqVKlGDx4MP369cPT05NZs2bh5ubG008/zeTJk5kzZ06uj8POzo7PP/+c4OBgzp8/z8svv8z48eOZPXs2TZo0YebMmUyePJnTp08bjwVg5MiRnDhxgiVLllCsWDFWrFjBk08+ydGjRylfvjx79+5lyJAhTJ8+ne7du7N27VqmTJmS6zhzSs5gQtwrJV6bO2jXl3ArQlvmVgQaDIN6QwrUjNLmSk038MaKo/x64AoAw5oH83qHytJnSOQ7SWl6qkxeZ5V9n3gnJMf/QDg4ODB//nyGDRvG3LlzqVOnDi1btqRPnz7UqFEj29e99tprhISEAPDqq6/St29fNm3aRNOmTQEYMmQI8+fPf6jjuLtGqXTp0rz33nu8+OKLzJ49GycnJ7y9vdHpdCZNepcuXWLevHlcunSJYsWKGWNdu3Yt8+bN4/3332fWrFk8+eSTjB8/HoAKFSqwa9cu1q5d+1DxPogkREJkSIzW5g/aMweSY7Rl3iWh6atQ+1lwtO3moNjENF786QC7z0dhp4O3u1XjuUZB1g5LiAKvV69edOrUib///ps9e/awZs0aZsyYwXfffZdtk9bdyZK/vz8A1atXN1kWERHxUHFt3LiR6dOnc+rUKeLi4khPTyc5OZnExMRs53w6evQoer2eChUqmCxPSUkxTpp58uRJevToYbK+cePGkhAJkecSImD3l/DP95CqtX9TuBw0GwM1ni4wF1t9GJeiEhk0fx/nIm/h7mTPl/3r0Lqin7XDEiLXXB3tOfFOiNX2bS4XFxfat29P+/btmTRpEkOHDmXKlCnZJkR3T1CY0Tx37zKDwWB8bmdnh1LKZBv3uwjuhQsX6Ny5My+99BLTpk3D19eXHTt2MGTIEFJTU7NNiBISErC3t+fAgQPY25u+DxlNatYiCZGwXTGXtf5BB3+E9GRtmX81aD4WqnQDO/NPWgXRwUs3GbZgP1G3UgnwcuGH5+tTpZiXtcMS4qHodLp83e+tSpUqFp17qGjRoly/ft34XK/Xc+zYMVq3bp1l+QMHDmAwGPjkk0+MV5r/5ZdfTMo4OTmh15v206pduzZ6vZ6IiAiaN2+e5bYrV67M3r17TZbt2bPH7GMyV/79NgiRG8mxcGYdnPhDu8+4tEbxetBiHFQIKZDXF8ut1UeuM+aXQ6SkG6hazIvvB9YnwNvF2mEJYTOioqJ46qmnGDx4MDVq1MDT05P9+/czY8YMunXrZrH9tGnThjFjxrB69WrKli3Lp59+SkxMTLbly5UrR1paGl988QVdunRh586dzJ0716RM6dKlSUhIYNOmTdSsWRM3NzcqVKhA//79GTBgAJ988gm1a9cmMjKSTZs2UaNGDTp16sQrr7xC06ZN+fjjj+nWrRvr1q3L8+YykGH3whYkRsO/P8Gip+GjcvDbMG3ovCENSjeHAX/A0I1Q8UlJhm67kZDCa8sOM2LxQVLSDbSt5McvwxtLMiTEI+bh4UHDhg357LPPaNGiBdWqVWPSpEkMGzaML7/80mL7GTx4MAMHDmTAgAG0bNmSMmXKZFs7BFCzZk0+/fRTPvzwQ6pVq8aiRYuYPn26SZkmTZrw4osv8swzz1C0aFFmzJgBwLx58xgwYABjx46lYsWKdO/enX/++YdSpUoB0KhRI7799ltmzZpFzZo1Wb9+PW+99ZbFjjU7OnVvo6HIJC4uDm9vb2JjY/HykqaCfCEhAk7+CSdXQujfoO6qti1SQWsSq9INAqpnvw0bpDcoFu+7xEdrTxGXnA7AkGbBvNFRRpKJ/Cs5OZnQ0FCCg4NxcZGkvqC53+drzu+3NJmJgkEpiDgJ57dqtT8Xd2GcSRq02aSrdNUusuon19nKyuHLMUz64xhHrsQCUCXQi3e7V6NuUCErRyaEEHlPEiKRf928AOe3Qeg2CN0OtyJN1xero9UCVe4ChctaJcT8ICYxlRnrTvPzvksoBZ7ODox9ogLPNgrCwV5a1YUQtkESIpF/JERoiU/oNi0Rirlout7BFYIaQ7n2WhLkU9I6ceYTBoPi1wNX+GDtKaJvpQLQo3ZxJnashJ+nNCsIIWyLJETi8RZ5Bo7+ol1MNeKE6To7B210WJmWENwSStQDh5xNa2/rTlyLY9Ifxzhw8SYAFfw9eKdbNRqVKWzlyIQQwjokIRKPn7jrcGy5lghdP2y6LqC6lvyUaQWlGoGzp1VCzK+OXY1lwa4LLD94BYMCNyd7RrUrz6CmwThK85gQwoZJQiQeD8mx2qiwI79ozWIZHaLtHKBcO6jWG8q2AXepwTBXarqBNceus2DXBQ5eijEu71QjkLc6VSbQ27YvSSKEECAJkbCm9BQ4u0GrCTq9FvQpd9aVbAQ1noIqPSQJyqXrsUks3nuJn/dd4kaC1kfI0V5Hh2qBDGxSWkaPCSHEXSQhEo9WQiT8t0GbJfrcZkiJu7OuaCWo/hRU7w2FSlstxPxMKcXu81Es3H2R9SfC0Ru0mjZ/L2f6NwyiT4OS0mFaCCGyIAmRyFtKQdgRLQE6sw6uHsBkfiDPYlC9F1R/WusfJDNF58q1mCQ2nAjnpz0XORuRYFzeMNiXgU1K076Kv/QREkKI+5CESFhe6i1tgsQz6+Dseoi/bro+sCaUD9GuG1asDtjJD7W50vQGDly8yZbTEWw9Fcnp8HjjOjcne3rULs6AxqWpGCCdzoUo6ObPn8+oUaPue+2xvPL8888TExNjsQvNTp06ld9//51Dhw5ZZHvmkIRI5J4+DaJD4cZpiDylDZGPPKXd9Kl3yjm6a6PCKoRA+SfAK9BqIedn4XHJbDsdyZbTEew4e4P4lHTjOjsd1C5ViM41AulVtwReLo5WjFQIYSnZJRxbt26ldevW3Lx5k2eeeYaOHTvmaHuWTp5mzZpFQbkCmCRE4sH06XDjjDYP0I0zd5KfqP/uXC3+Xj5BUOFJqPAEBDUDR+m3Yq6ohBSOXYtjX2gUW05FcuJ6nMl6X3cnWlUoSqtKfrQoXwQfNycrRSqEsCZXV1dcXR/taFG9Xo9Op8Pb2/uR7jcvSUIkTKWnQuRJbf6fjFvYMUhPyrq8ozsUKa91iC5aQbv3qwyFgqU/UA4ppYiIT+HY1ViOXY3j6NVYjl+L5Xpsskk5nQ5qlPChdcWitKroR43i3tjJBVeFsHn31vocPnyYUaNGsX//fnQ6HeXLl+frr78mISGBQYMGAaC7fX6eMmUKU6dO5ebNm7z66qv8+eefpKSk0LJlSz7//HPKly9vso8ff/yR119/nTNnzvDff/8xdepUkxosg8HAxx9/zDfffMPly5fx9/dn+PDhvPnmmwBMmDCBFStWcOXKFQICAujfvz+TJ0/G0dH6tdqSENmy9BQIP6YlPdcOafcRJ0ybuzI4eWqJTtGKt5Of2wmQVwnpA2SGuOQ0LkUlcik6kZPX4zh2NZajV+O4kZCSZfkyRdypUcKbFhWK0qJCUYp4yEzcQliEUpCWaJ19O7rl6T+M/fv3p3bt2syZMwd7e3sOHTqEo6MjTZo0YebMmUyePJnTp08D4OHhAWhNc2fPnmXlypV4eXkxYcIEOnbsyIkTJ4zJSmJiIh9++CHfffcdhQsXxs/PL9O+J06cyLfffstnn31Gs2bNuH79OqdOnTKu9/T0ZP78+RQrVoyjR48ybNgwPD09GT9+fJ69HzklCZGtUEq79teV/bdv/2ijv7JKfly8tY7PgbXu3PuWkcQnB9L1Bq7HJnM5Wkt6LkUncjE60fg8JjHrJkY7HZTz86BaMW+qFddulQM98ZS+QELkjbREeL+Ydfb9xjVwcs9x8VWrVhkTlwx6vT7b8pcuXWLcuHFUqlQJwFjLA+Dt7Y1OpyMgIMC4LCMR2rlzJ02aNAFg0aJFlCxZkt9//52nnnoKgLS0NGbPnk3NmjWz3G98fDyzZs3iyy+/ZODAgQCULVuWZs2aGcu89dZbxselS5fmtddeY8mSJZIQPWpfffUVH330EWFhYdSsWZMvvviCBg0aWDusvJEcC1cPwtX9d5KgxBuZy7n6QrHatxOfmlCsltb/R5q7jPQGxc3EVG4kpHAjPpWoWylExqdwI+H2soxbvPY83XD/DoaF3Z0o4etGRX8PqhX3pmoxb6oEeuHqZP+IjkgIkZ+0bt2aOXPmmCzbu3cvzz77bJblx4wZw9ChQ1m4cCHt2rXjqaeeomzZstlu/+TJkzg4ONCwYUPjssKFC1OxYkVOnjxpXObk5ESNGjXuu52UlBTatm2bbZmlS5fy+eefc+7cORISEkhPT8fLyyvb8o+SzSRES5cuZcyYMcydO5eGDRsyc+ZMQkJCOH36dJbVfo81gx5uRULcVYi7pl37y/j4GsRdgZsXMZnvB8DOEQJraBdELVEfStQtkH19lFKk6RVJqXqS0vQkpqaTlKYnOU1PYqreuPxWip745DTik9OJy7hPyvw8ITUdcwZRONnbUcLXlVK+bsZbybvuPZxt5s9OiMeXo5tWU2OtfZvB3d2dcuXKmSy7cuVKtuWnTp1Kv379WL16NWvWrGHKlCksWbKEHj165CrcDK6ursa+R9mtv5/du3fTv39/3n77bUJCQvD29mbJkiV88sknDxWXpdjMmfnTTz9l2LBhxg5lc+fOZfXq1fzwww+8/vrrVolJn5ZC5IXj2KUloEuJR5eagF1qPLrUeHQpCdil3Xlul5qALjkWh1vXsU+MQGdIf+D2Uz1Lkexfi0S/OiQVrU1S4coY7LXRXkqBSlSoW7EowKCUtkwp7blBuzdZphQGlVFWYTCYLst4rJRCb1CkGxQGg0Kvbt8bFHoFeoMB/e3XpusV6QYD6QZFut5Amj7jtXcep+kNpOsVqXoDqenaLcX4WG+yPDXdQHK6wThDsyUVcnOkiIezdvN0poiHE0U8nCnq4UwRz9uPPZ3x93SRzs5CPO50OrOarfKbChUqUKFCBUaPHk3fvn2ZN28ePXr0wMnJKVNzW+XKlUlPT2fv3r3GJrOoqChOnz5NlSpVcrzP8uXL4+rqyqZNmxg6dGim9bt27SIoKMjYwRrg4sWLuTxCy7OJhCg1NZUDBw4wceJE4zI7OzvatWvH7t27M5VPSUkhJeVOJ9e4uLhMZSzhZsRlAha1ztVr9UpHOIUIV75cV76EKV/CVCHCVGGuK19CVSA3kr0hMuMVt4D9lgo937C30+HmaI+r0+2b4517Nyd7vFwc8XRxwMv19r2LI54ujni5Omj3Ltq9j5ujzPQshHjsJSUlMW7cOHr37k1wcDBXrlzhn3/+oVevXoDWbychIYFNmzZRs2ZN3NzcKF++PN26dWPYsGF8/fXXeHp68vrrr1O8eHG6deuW4327uLgwYcIExo8fj5OTE02bNiUyMpLjx48zZMgQypcvz6VLl1iyZAn169dn9erVrFixIq/eCrPZREJ048YN9Ho9/v7+Jsv9/f1Ner9nmD59Om+//XbeB+bkRbTyJAFXEnDjFq4k4Hr73o1EXEhQt9fp3EjAjUgKEU5honTeGLAH3Z3hkxk1mTq0ZRmXRNWW69DpMtaBzuS59thOd2eZ3e0Hdjqd8TV2Op12s8soq8POuFzbjv3tbdjb6bSbTofd7Xt7O+2xg522HXs7sLezw9Feh4OdHQ722joHezsc7XTY2+twzFhub4ezvR1ODrdvdz+++7m9Ha5O9rjcTngkiRFC2BJ7e3uioqIYMGAA4eHhFClShJ49exp/05o0acKLL77IM888Q1RUlHHY/bx583j11Vfp3LkzqamptGjRgr/++svs4fCTJk3CwcGByZMnc+3aNQIDA3nxxRcB6Nq1K6NHj2bkyJGkpKTQqVMnJk2axNSpUy39NuSKThWUKSbv49q1axQvXpxdu3bRuHFj4/Lx48ezbds29u7da1I+qxqikiVLEhsb+9h0/hJCCJEzycnJhIaGEhwcjIuLTBJb0Nzv842Li8Pb2ztHv982UUNUpEgR7O3tCQ8PN1keHh5uMvQwg7OzM87OMt+LEEIIYStsoj3BycmJunXrsmnTJuMyg8HApk2bTGqMhBBCCGGbbKKGCLR5GQYOHEi9evVo0KABM2fO5NatW8ZRZ0IIIYSwXTaTED3zzDNERkYyefJkwsLCqFWrFmvXrs3U0VoIIYQQtsdmEiKAkSNHMnLkSGuHIYQQQojHjE30IRJCCCFsYFC1TbLU5yoJkRBCiALt7qu1i4InNVW7SLm9/cNdD9KmmsyEEELYHnt7e3x8fIiIiADAzc3tvtfkEvmHwWAgMjISNzc3HBweLqWRhEgIIUSBlzHnXEZSJAoOOzs7SpUq9dBJriREQgghCjydTkdgYCB+fn6kpaVZOxxhQU5OTtjZPXwPIEmIhBBC2Ax7e/uH7msiCibpVC2EEEIImycJkRBCCCFsniREQgghhLB50ocoBzImfYqLi7NyJEIIIYTIqYzf7ZxM3igJUQ7Ex8cDULJkSStHIoQQQghzxcfH4+3tfd8yOiVzmT+QwWDg2rVreHp6Wnwyr7i4OEqWLMnly5fx8vKy6LaFZchn9PiTz+jxJ5/R462gfj5KKeLj4ylWrNgDh+ZLDVEO2NnZUaJEiTzdh5eXV4H6EhZE8hk9/uQzevzJZ/R4K4ifz4NqhjJIp2ohhBBC2DxJiIQQQghh8yQhsjJnZ2emTJmCs7OztUMR2ZDP6PEnn9HjTz6jx5t8PtKpWgghhBBCaoiEEEIIISQhEkIIIYTNk4RICCGEEDZPEiIhhBBC2DxJiKzoq6++onTp0ri4uNCwYUP27dtn7ZBs2vbt2+nSpQvFihVDp9Px+++/m6xXSjF58mQCAwNxdXWlXbt2nD171jrB2qDp06dTv359PD098fPzo3v37pw+fdqkTHJyMiNGjKBw4cJ4eHjQq1cvwsPDrRSx7ZkzZw41atQwTu7XuHFj1qxZY1wvn8/j5YMPPkCn0zFq1CjjMlv+jCQhspKlS5cyZswYpkyZwsGDB6lZsyYhISFERERYOzSbdevWLWrWrMlXX32V5foZM2bw+eefM3fuXPbu3Yu7uzshISEkJyc/4kht07Zt/2/v/mOirv84gD/vDg6hOw8QdhcxfpRn6BRFGIEuU2Fs1Cii1XLMDvzRpIN5OIoxOddmy1ptgTmbowX+A2dkzMVmRkZsBRQ/ugIxE6ZS84BZQoDl0d37+4fjsy6zfTHgc+3zfGzvjc/787n7PPW1g9c+7/vctcFqtaKzsxMtLS2YmZlBVlYWpqenpWNKS0vx0UcfobGxEW1tbbh69Sry8vJkTK0s0dHReO2119DT04Pu7m5s3boVTzzxBM6dOweA9fEnXV1dOHbsGBITE33mFV0jQbJITU0VVqtV2vZ4PCIqKkocOnRIxlQ0C4BoamqStr1erzCZTOKNN96Q5sbHx0VQUJBoaGiQISGNjY0JAKKtrU0IcasegYGBorGxUTrm/PnzAoDo6OiQK6bihYWFiXfffZf18SOTk5PCbDaLlpYW8cgjj4i9e/cKIfga4hUiGbjdbvT09CAzM1OaU6vVyMzMREdHh4zJ6E4uXbqEkZERn5oZDAY89NBDrJlMJiYmAADh4eEAgJ6eHszMzPjUKCEhATExMayRDDweDxwOB6anp5Gens76+BGr1YrHHnvMpxYAX0P8clcZXLt2DR6PB0aj0WfeaDTi+++/lykV/ZORkREA+Nuaze6jxeP1emGz2bBx40asXr0awK0aabVahIaG+hzLGi2uvr4+pKen4/fff4dOp0NTUxNWrVoFp9PJ+vgBh8OB3t5edHV13bZP6a8hNkRE9J9jtVrR39+PL774Qu4o9BcPPvggnE4nJiYm8MEHH8BisaCtrU3uWATgxx9/xN69e9HS0oIlS5bIHcfvcMlMBhEREdBoNLe9c390dBQmk0mmVPRPZuvCmsmvuLgYzc3NaG1tRXR0tDRvMpngdrsxPj7uczxrtLi0Wi2WL1+O5ORkHDp0CGvXrkV1dTXr4wd6enowNjaG9evXIyAgAAEBAWhra8Phw4cREBAAo9Go6BqxIZKBVqtFcnIyzp49K815vV6cPXsW6enpMiajO4mPj4fJZPKp2a+//oqvvvqKNVskQggUFxejqakJn332GeLj4332JycnIzAw0KdGFy5cwPDwMGskI6/Xi5s3b7I+fiAjIwN9fX1wOp3SSElJQX5+vvSzkmvEJTOZ7Nu3DxaLBSkpKUhNTUVVVRWmp6dRWFgodzTFmpqawuDgoLR96dIlOJ1OhIeHIyYmBjabDa+88grMZjPi4+Nht9sRFRWF3Nxc+UIriNVqRX19PU6dOgW9Xi+9p8FgMCA4OBgGgwE7d+7Evn37EB4ejqVLl6KkpATp6elIS0uTOb0yVFRUIDs7GzExMZicnER9fT0+//xznDlzhvXxA3q9XnrP3ax77rkHy5Ytk+YVXSO5b3NTsrffflvExMQIrVYrUlNTRWdnp9yRFK21tVUAuG1YLBYhxK1b7+12uzAajSIoKEhkZGSICxcuyBtaQf6uNgBEbW2tdMxvv/0mXnjhBREWFiZCQkLEk08+KVwul3yhFWbHjh0iNjZWaLVaERkZKTIyMsQnn3wi7Wd9/M+fb7sXQtk1UgkhhEy9GBEREZFf4HuIiIiISPHYEBEREZHisSEiIiIixWNDRERERIrHhoiIiIgUjw0RERERKR4bIiIiIlI8NkREJKuCggJZPu27rq4OKpUKKpUKNptt0c///5jN99dvHyei+ceGiIgWzOwf9DuNl19+GdXV1airq5Ml39KlS+FyuXDw4EFpbvPmzVK+oKAg3HfffcjJycGHH3646PlcLheqqqoW/bxESsSGiIgWjMvlkkZVVZXUgMyOsrIyGAwG2a6AqFQqmEwm6PV6n/ndu3fD5XJhaGgIJ0+exKpVq/Dss8/i+eefX9R8JpMJBoNhUc9JpFRsiIhowZhMJmkYDAapAZkdOp3utiWzzZs3o6SkBDabDWFhYTAajaipqZG+/Fiv12P58uU4ffq0z7n6+/uRnZ0NnU4Ho9GI7du349q1a3eVOyQkBCaTCdHR0UhLS8Prr7+OY8eOoaamBp9++ql0XHl5OVasWIGQkBDcf//9sNvtmJmZAQBcvnwZarUa3d3dPs9dVVWF2NhYeL1eXL9+Hfn5+YiMjERwcDDMZjNqa2vvKjMR/TtsiIjI7xw/fhwRERH4+uuvUVJSgqKiIjz99NPYsGEDent7kZWVhe3bt+PGjRsAgPHxcWzduhVJSUno7u7Gxx9/jNHRUTzzzDPzlslisSAsLMxn6Uyv16Ourg4DAwOorq5GTU0N3nrrLQBAXFwcMjMzb2twamtrUVBQALVaDbvdjoGBAZw+fRrnz5/HO++8g4iIiHnLTET/PzZEROR31q5di8rKSpjNZlRUVGDJkiWIiIjA7t27YTabceDAAfz888/47rvvAABHjhxBUlISXn31VSQkJCApKQnvvfceWltb8cMPP8xLJrVajRUrVuDy5cvSXGVlJTZs2IC4uDjk5OSgrKwM77//vrR/165daGhowM2bNwEAvb296OvrQ2FhIQBgeHgYSUlJSElJkRqonJyceclLRHPDhoiI/E5iYqL0s0ajwbJly7BmzRppzmg0AgDGxsYAAN9++y1aW1uh0+mkkZCQAAAYGhqat1xCCKhUKmn7xIkT2Lhxo7T8V1lZieHhYWl/bm4uNBoNmpqaANy6s23Lli2Ii4sDABQVFcHhcGDdunV46aWX0N7ePm9ZiWhu2BARkd8JDAz02VapVD5zs02J1+sFAExNTSEnJwdOp9NnXLx4EZs2bZqXTB6PBxcvXkR8fDwAoKOjA/n5+Xj00UfR3NyMb775Bvv374fb7ZYeo9Vq8dxzz6G2thZutxv19fXYsWOHtD87OxtXrlxBaWkprl69ioyMDJSVlc1LXiKamwC5AxAR/Vvr16/HyZMnERcXh4CAhfm1dvz4cVy/fh1PPfUUAKC9vR2xsbHYv3+/dMyVK1due9yuXbuwevVqHD16FH/88Qfy8vJ89kdGRsJiscBiseDhhx/Giy++iDfffHNB/g1EdGe8QkRE/3lWqxW//PILtm3bhq6uLgwNDeHMmTMoLCyEx+OZ8/PduHEDIyMj+Omnn9DZ2Yny8nLs2bMHRUVF2LJlCwDAbDZjeHgYDocDQ0NDOHz4sLQ09mcrV65EWloaysvLsW3bNgQHB0v7Dhw4gFOnTmFwcBDnzp1Dc3MzVq5ceff/EUR019gQEdF/XlRUFL788kt4PB5kZWVhzZo1sNlsCA0NhVo9919zNTU1uPfee/HAAw8gLy8PAwMDOHHiBI4ePSod8/jjj6O0tBTFxcVYt24d2tvbYbfb//b5du7cCbfb7bNcBtxaUquoqEBiYiI2bdoEjUYDh8Mx57xE9O+phBBC7hBERIutrq4ONpsN4+PjC36ugwcPorGxUborbi4WMyeRkvEKEREp1sTEBHQ6HcrLyxfk+aemptDf348jR46gpKRkzo/X6XTYs2fPAiQjor/iFSIiUqTJyUmMjo4CAEJDQxfkAxELCgrQ0NCA3Nxc1NfXQ6PRzOnxg4ODAG599MDs3W1EtDDYEBEREZHiccmMiIiIFI8NERERESkeGyIiIiJSPDZEREREpHhsiIiIiEjx2BARERGR4rEhIiIiIsVjQ0RERESKx4aIiIiIFO9/7QnNOPwJ1/YAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "population, contact_frequency = res.x\n", "result = model.run(params={'total_population':population,\n", " 'contact_frequency':contact_frequency},\n", " return_columns=['population_infected_with_ebola'],\n", " return_timestamps=list(data.index.values))\n", "\n", "plt.plot(result.index, result['population_infected_with_ebola'], label='Simulated')\n", "plt.plot(data.index, data['Cumulative Cases'], label='Historical');\n", "plt.xlabel('Time [Days]')\n", "plt.ylabel('Cumulative Infections')\n", "plt.title('Model fit to Sierra Leone Ebola historical infections data')\n", "plt.legend(loc='lower right')\n", "plt.text(2,9000, 'RMSE: 7.5% of Max', color='r', fontsize=12)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " fun: 22200248.07232056\n", " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([2.60749363e+00, 3.06777631e+03])\n", " message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 93\n", " nit: 10\n", " njev: 31\n", " status: 0\n", " success: True\n", " x: array([8.82124175e+03, 8.20469365e+00])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.07505469143880455" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sqrt(res.fun/len(data))/data['Cumulative Cases'].max()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.10 64-bit", "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.8.10" }, "vscode": { "interpreter": { "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1" } } }, "nbformat": 4, "nbformat_minor": 1 }