Fudan_CCDC_Japan_Optuna.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fudan CCDC model (time delay dynamical system):\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\mathrm{d} I}{\\mathrm{d} t} &= r I_{0}(t) \\\\\n",
    "\\frac{\\mathrm{d} J}{\\mathrm{d} t} &= r \\int_{-\\infty}^{t} f_{4}(t-s) I_{0}(s) \\mathrm{d} s \\\\\n",
    "\\frac{\\mathrm{d} G}{\\mathrm{d} t} &= \\ell(t) \\left[ \\int_{-\\infty}^{t} f_{2}(t-s) I_{0}(s) d s - \\int_{-\\infty}^{t} f_{4}(t-s) I_{0}(s) \\mathrm{d} s \\right] \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "- $I(t)$: the cumulative number of the infected at time $t$;\n",
    "\n",
    "- $J(t)$: the cumulative number of confirmed in hospitalization stage at time $t$;\n",
    "\n",
    "- $G(t)$: the number of infected and isolated but yet diagnosed cases at time $t$. They are infected in fact, but are not conﬁrmed by the hospital, then not appeared in the infected list of CCDC;\n",
    "\n",
    "- $I_{0}(t) := I(t) - J(t) - G(t)$: the number of unconfirmed/unobserved infected at time $t$;\n",
    "\n",
    "\n",
    "<div align=\"center\">\n",
    "Fudan CCDC's graphical model<br>\n",
    "<img src=\"fudanCCDC.png\" alt=\"png\" title=\"Fudan CCDC\" width=20%>\n",
    "</div>\n",
    "\n",
    "References:\n",
    "\n",
    "[The reproductive number R0 of COVID-19 based on estimate of a statistical time delay dynamical system][1]\n",
    "\n",
    "[CoVID-19 in Japan: What could happen in the future?][2]\n",
    "\n",
    "[Qiita記事: 中国の最新論文の方式で日本のコロナウィルスの感染者数を予測してみた][3]\n",
    "\n",
    "Data sources:\n",
    "\n",
    "[厚生労働省: 新型コロナウイルス感染症の現在の状況について（３月８日12時時点版）][4]\n",
    "\n",
    "\n",
    "[1]:https://www.medrxiv.org/content/10.1101/2020.02.17.20023747v2\n",
    "[2]:https://www.medrxiv.org/content/10.1101/2020.02.21.20026070v2\n",
    "[3]:https://qiita.com/KEN_KEN2/items/a1aae5eb4129e5fbd91f\n",
    "[4]:https://www.mhlw.go.jp/stf/newpage_10032.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzgAAAG5CAYAAABC/5HWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXRU9f3/8dcnkx0SCGHf17AT1ICIdfu6i7iD1tqiGBGXWrV1aetWtBaXqnUvLu2vVqsggqAoYkEFlwpWkrAOIUAIZAIJW4Akk8x8fn9koAERJpDkzvJ8nMMxM/femddwjklefO59X2OtFQAAAABEghinAwAAAABAQ6HgAAAAAIgYFBwAAAAAEYOCAwAAACBiUHAAAAAARAwKDgAAAICIQcEBADQKY8zLxpj7nc4BAIguhvvgAED4MMasl9ROUo0kn6QVkv4haYq11u9gtAZT5zP6JO2R9JGkW621u53MdSTGmL9LKrLW3ud0FgCIZqzgAED4GW2tTZHUTdJkSfdIes3ZSA1utLW2uaTjJWVJqndpMMbENniqRhRueQEgVFFwACBMWWt3WmtnSbpS0jhjzCBjzDBjTIkxxrVvP2PMZcaYnMDXDxljphpj/mGMKTfGLDfGZNXZ915jzNrAthXGmEvrbLvWGPOlMeZpY8wOY0yBMWZk4PmNxpgtxphxdfb/uzHmkTqPLzbGLDXG7Aq8x3lBfMZNql3BGRR4jeuMMSsD+QqMMTfWef3TjTFFxph7jDEeSX8zxqQZYz4wxmw1xmwPfN25zjGfGWMeMcZ8ZYzZbYyZbYxJN8a8Gci52BjTvc7+/Ywx84wx24wxq40xYwPPT5D0M0l373udwPMdjTHTA++/zhhzW53XesgY864x5p/GmF2SrjXGDDfGLAm8d4kx5qkj/R0BAA5EwQGAMGet/VZSkaRTrLWLJZVJOqfOLj9X7Wls+1wk6W1JLSXNkvR8nW1rJZ0iqYWkP0j6pzGmQ53tJ0rKlZQu6a3A6wyT1FvSNZKeN8Y0PzijMWZ4IMNdgfc9VdL6I302Y0wXSRdI+j7w1BZJF0pKlXSdpKeNMcfXOaS9pFaqXd2aoNqfc38LPO4qqeKgzytJV6n276iTpF6Svg4c00rSSkkPBrI0kzQv8LnbBo570RgzwFo7RdKbkh631ja31o42xsRImi0pJ/DaZ0q63Rhzbp33vljSu4G/kzcl/UXSX6y1qYEsU4/0dwQAOBAFBwAiw2bV/kIuSf9PtWVDxphWks5V7S/l+yyy1s6x1vokvSEpc98Ga+00a+1ma63fWvuOpDWShtc5dp219m+BY9+R1EXSJGttlbX2E0le1Zadg10v6XVr7bzAa2+y1q46zOeZaYzZIWmRpM8lPRrI96G1dq2t9bmkT1RbyPbxS3owkKfCWltmrZ1urd1rrS2X9EdJpx30Xn8LvOZO1a4WrbXWfmqtrZE0TdJxgf0ulLQ+8PlrrLXfS5ouacyPfIZhktpYaydZa73W2gJJr6i2GO3ztbV2ZuDvpEJStaTexpjW1trd1tpvDvN3BAA4BM73BYDI0EnStsDX/5S0MrDiMFbSQmttcZ19PXW+3isp0RgTa62tMcb8QtKdkroHtjeX1LrO/iV1vq6QJGvtwc/9YAVHtUVoTj0+zyXW2k8PftIYc75qV1QyVPuPdMmS8ursstVaW1ln/2RJT0s6T1Ja4OkUY4wrUNIO9Zl+7PN0k3RioHjtE6vakngo3SR1PGh/l6SFdR5vPOiY6yVNkrTKGLNO0h+stR/8yOsDAA6BggMAYc4YM0y1BWeRVHvdijHma0mXqfbUq5eCfJ1uql1hOFO1Kws+Y8xSSaYBYm5U7SlXR80Yk6DaFZNfSHrfWlttjJl5UL6DR4P+WlJfSSdaaz3GmKGqPd3taD7TRkmfW2vP/pHtB7/3RtWuePU5zGsecIy1do2knwZOb7tM0rvGmHRr7Z6jyAsAUYlT1AAgTBljUo0xF6r2Oph/WmvrrmT8Q9LdkgZLei/Il2ym2l+4twZe/zoFLu5vAK9Jus4Yc6YxJsYY08kY06+erxEvKSGQryawmnPO4Q9RimpXYXYETtd7sL7B6/hAUoYx5ufGmLjAn2HGmP6B7SWSetbZ/1tJ5YGhB0nGGJcJDIL4sTcwxlxjjGkTGPm9b+UnIsZ/A0BToeAAQPiZbYwpV+0Kwe8lPaXaC+7rmqHaU6RmWGv3BvOi1toVkv6s2ovsS1Rbjr5siMCBQQjXqfZ0sZ2qva6mWz1fo1zSbaq98H67pKtVOyThcJ6RlCSpVNI3kj6uV/Afvv85qr2GZrNqT/V7TLWlS6otcQMCE+ZmBk6Bu1DSUEnrAhleVe0Ahx9znqTlxpjdqh04cFXg2hwAQJC40ScARChjzFpJNx7qWhYAACIVKzgAEIGMMZer9nSz+U5nAQCgKTFkAAAijDHmM0kDJP08cC0HAABRg1PUAAAAAEQMTlEDAAAAEDFC7hS11q1b2+7duzsdAwAAAEAI++6770qttW0Ofj7kCk737t21ZMkSp2MAAAAACGHGmA2Hep5T1AAAAABEDAoOAAAAgIhBwQEAAAAQMULuGpxDqa6uVlFRkSorK52O0mgSExPVuXNnxcXFOR0FAAAACFthUXCKioqUkpKi7t27yxjjdJwGZ61VWVmZioqK1KNHD6fjAAAAAGErLE5Rq6ysVHp6ekSWG0kyxig9PT2iV6gAAACAphAWBUdSxJabfSL98wEAAABNIWwKDgAAAAAcCQUHAAAAQMSg4NTDs88+q/79+8sYoyFDhmjw4MEaOXKkcnJy9u9TUVGh0047TT6fT+vXr9dbb721f1teXp6uvfZaB5IDAAAA0YGCUw8vvvii5s2bpy+//FKff/658vLydP/992vChAn793n99dd12WWXyeVy/aDgDB48WEVFRSosLHQiPgAAABDxwmJMdF23z1ympZt3NehrDu2YqmcuGXTYfSZOnKiCggKdf/75Gj9+vEaOHClJGjFihIqKivbv9+abb+4vNffee69WrlypoUOHaty4cbrjjjs0evRovf3227r77rsb9DMAAAAACHIFxxhznjFmtTEm3xhz7yG232mMWWGMyTXG/NsY063ONp8xZmngz6yGDN+UXn75ZXXs2FELFizQHXfcsf/51157Teeff74kyev1qqCgQN27d5ckTZ48WaeccoqWLl26/5isrCwtXLiwyfMDAAAA0eCIKzjGGJekFySdLalI0mJjzCxr7Yo6u30vKctau9cYc5OkxyVdGdhWYa0d2lCBj7TS0pQWLFig1157TYsWLZIklZaWqmXLloc9pm3bttq8eXNTxAMAAACiTjCnqA2XlG+tLZAkY8zbki6WtL/gWGsX1Nn/G0nXNGTIUJSbm6vs7Gx99NFHSk9PlyQlJSUd8WadlZWVSkpKaoqIAAA0OVvjVZXH7XQMAI3AxMQqoWM/p2McUTAFp5OkjXUeF0k68TD7Xy/pozqPE40xSyTVSJpsrZ1Z75QhprCwUJdddpneeOMNZWRk7H8+LS1NPp9PlZWVSkxMVEpKisrLyw841u12a9Cg0FmFAgCgIXnevF3b57/kdAwAjSA2rZMynik68o4Oa9AhA8aYayRlSTqtztPdrLWbjDE9Jc03xuRZa9cedNwESRMkqWvXrg0ZqVFMmjRJZWVluvnmmyVJsbGxWrJkiSTpnHPO0aJFi3TWWWdpyJAhcrlcyszM1LXXXqs77rhDCxYs0KhRo5yMDwBAo6nc8L0SumaqzUX3OR0FQAMzceFxFlIwBWeTpC51HncOPHcAY8xZkn4v6TRrbdW+5621mwL/LTDGfCbpOEkHFBxr7RRJUyQpKyvL1u8jNJ3169dLkl599VW9+uqrh9znlltu0dNPP62zzjpLcXFxmj9//v5tVVVVWrJkiZ555pmmiAsAQJPzetxKHT5GqcOucDoKgCgVzBS1xZL6GGN6GGPiJV0l6YBpaMaY4yT9VdJF1totdZ5PM8YkBL5uLelk1bl2JxIdf/zxOuOMM+Tz+X6wrbCwUJMnT1ZsbNhN5wYA4IhqdpfJt2eb4ttnHHlnAGgkR/xN21pbY4y5VdJcSS5Jr1trlxtjJklaYq2dJekJSc0lTTPGSFKhtfYiSf0l/dUY41dtmZp80PS1iDR+/PhDPt+nTx/16dOnidMAANA0vIHhAhQcAE4KainBWjtH0pyDnnugztdn/chxX0kafCwBAQBAeKDgAAgFQd3oEwAA4Ei8HrfkilV86x5ORwEQxSg4AACgQVR53Ipv01MmNs7pKACiGAUHAAA0CG/xak5PA+A4Ck49PPvss+rfv79+9rOfSZIWL16s2NhYvfvuu/v3KS4u1oUXXihJWrp0qebM+d+lSx988IEeeOABAQAQaazfL2/JGgoOAMdRcOrhxRdf1Lx58/Tmm2/K5/Ppnnvu0TnnnHPAPk899ZRuuOEGST8sOKNGjdLs2bO1d+/eJs0NAEBjq9leJFtdqQQKDgCHhd0NWW7/z/tauu0H9xk9JkNbddIzJ1582H0mTpyogoICnX/++Ro/fryMMbr88su1ePHiA/abPn26HnnkEXm9Xj3wwAOqqKjQokWL9Nvf/lZXXnmlTj/9dH3wwQcaO3Zsg34GAACcVMUENQAhghWcIL388svq2LGjFixYoLFjx2rGjBm66aabDthn3bp1SktLU0JCguLj4zVp0iRdeeWVWrp0qa688kpJUlZWlhYuXOjERwAAoNHsHxHdjoIDwFlht4JzpJWWpnD77bfrscceU0zMgf2wuLhYbdq0Oeyxbdu21ebNmxszHgAATc7rccvEJys2raPTUQBEubArOKFgyZIluuqqqyRJpaWlmjNnjmJjY9WtWzdVVlYe9tjKykolJSU1RUwAAJqM1+NWQvsMGWOcjgIgylFwjsK6dev2f33ttdfqwgsv1CWXXKI9e/Zo/fr1+7elpKSovLz8gGPdbrcGDRrUVFEBAGgSVR63krqf4HQMAOAanIbUrFkz9erVS/n5+ZKkM844QytWrNDQoUP1zjvvSJIWLFigUaNGORkTAIAGZWu8qt66jgEDAEICKzj1UHd1Zp+///3vBzy+9dZb9fe//12PPPKIWrVqdcCUtZKSElVUVGjw4MGNnBQAgKbj3VIgWT8FB0BIoOA0sEsvvVRlZWWH3FZYWKg///nPTZwIAIDGtW+CWkL7vg4nAQAKTqPIzs4+5PPDhg1r4iQAADS+/90Dp4/DSQCAa3AAAMAx8npWy5XSRq5maU5HAQAKDgAAODZej5vrbwCEDAoOAAA4JvvugQMAoYCCAwAAjpqvYpdqdnpYwQEQMig4QWrevLmk2lHRb7311gHbvv/+e11//fWSpM8++0xfffXV/m3PP/+8Xn/99aYLCgBAE/J61kgSBQdAyKDg1NOhCs6jjz6q2267TdIPC8748eP13HPPNWlGAACainf/BDUKDoDQEHZjoj1v3q7KwqUN+pqJXYeq/c+eCWrfe++9VytXrtTQoUM1btw4ZWdnKzc3V5mZmVq/fr1efvlluVwu/fOf/9Rzzz2nU045Rd27d9e3336r4cOHN2huAACc5vW4JWMU37aX01EAQFIYFhynTZ48WU8++aQ++OADSdKCBQs0aNAgSVL37t01ceJENW/eXL/5zW/2H5OVlaWFCxdScAAAEaeqxK249K6KiU9yOgoASArDghPsSktTKS4uVps2bQ67T9u2bbVq1aomSgQAQNPxetyKb8fpaQBCB9fgHKOkpCRVVlYedp/KykolJfEvWwCAyGKt5R44AEIOBaeeUlJSVF5evv9x//79lZ+f/6PbJcntdu8/jQ0AgEjh27VF/opdSujQ1+koALAfBaeehgwZIpfLpczMTD399NPq16+fdu7cub/UjB49WjNmzNDQoUO1cOFCSdKXX36ps88+28nYAAA0uComqAEIQWF3DY5Tdu/eLUmKi4vT/PnzD9g2fvx4vfPOO8rOzlZGRoZyc3P3b/v+++81cOBApaenN2leAAAam7d4tSQKDoDQwgpOA7jpppuUkJBwyG2lpaV6+OGHmzgRAACNz1vilomNV1x6V6ejAMB+YbOCY62VMcbpGIeUmJion//854fcFuypadbahowEAECj83rcim/bWybG5XQUANgvLFZwEhMTVVZWFrElwFqrsrIyJSYmOh0FAICgVTFBDUAICosVnM6dO6uoqEhbt251OkqjSUxMVOfOnZ2OAQBAUKzfJ29JvlKGjnY6CgAcICwKTlxcnHr06OF0DAAAEFBdukHyVbOCAyDkhMUpagAAILR4GRENIERRcAAAQL3tuwdOAgUHQIih4AAAgHrzetyKSUqVK7Wt01EA4AAUHAAAUG9ej1vx7fuG7C0cAEQvCg4AAKg3b4mb09MAhCQKDgAAqBe/t0LVZYUMGAAQkig4AACgXrwl+ZK1FBwAIYmCAwAA6oUR0QBCGQUHAADUy/6C066Pw0kA4IcoOAAAoF6qPG7FtuwgV1KK01EA4AcoOAAAoF5qR0RzehqA0ETBAQAA9eItcSu+HQUHQGii4AAAgKD5dm+Tr7yUFRwAIYuCAwAAglZVskaSuMkngJBFwQEAAEHbP0GtQ1+HkwDAoVFwAABA0LwetxTjUnybHk5HAYBDouAAAICgeT1uxbfpIRMb73QUADgkCg4AAAgaI6IBhDoKDgAACIr1+1XlYUQ0gNBGwQEAAEGp2bFZ1ruXFRwAIY2CAwAAgrJ/ghoFB0AIo+AAAICgVAUKDvfAARDKKDgAACAoXo9bJj5JsWmdnI4CAD+KggMAAILi9bgV366PTAy/PgAIXXyHAgAAQWFENIBwQMEBAABHZGuq5d1awPU3AEIeBQcAAByRt3Sd5Pcpvn1fp6MAwGFRcAAAwBExIhpAuKDgAACAI/IyIhpAmAiq4BhjzjPGrDbG5Btj7j3E9juNMSuMMbnGmH8bY7rV2TbOGLMm8GdcQ4YHAABNw+txy9U8Xa7mrZyOAgCHdcSCY4xxSXpB0vmSBkj6qTFmwEG7fS8py1o7RNK7kh4PHNtK0oOSTpQ0XNKDxpi0hosPAACaQlXxak5PAxAWYoPYZ7ikfGttgSQZY96WdLGkFft2sNYuqLP/N5KuCXx9rqR51tptgWPnSTpP0r+OPToAADhYSXmVbp+5TO/mFstvbYO97qdlufombqju+83sBntNAOGlU4tEFd5/ttMxjiiYgtNJ0sY6j4tUuyLzY66X9NFhjv3B7Y+NMRMkTZCkrl27BhEJAADUZa3V3xdv1K9nrdAer08TRnRTq+S4Bnnt2Oo9ajd9m7pnDNHvBvZpkNcEEH5SE4OpDs5r0JTGmGskZUk6rT7HWWunSJoiSVlZWQ33z00AAESBtaV7dOO7ufr3mlL9pEcrvTJmiPq1S2mw16/Y8L3WTZdGnXayfjqsX4O9LgA0hmAKziZJXeo87hx47gDGmLMk/V7SadbaqjrHnn7QsZ8dTVAAAHCgGp9fT39RoAfnrlZsTIxeunywJozoppgY06Dvs39EdDuuwQEQ+oIpOIsl9THG9FBtYblK0tV1dzDGHCfpr5LOs9ZuqbNprqRH6wwWOEfSb485NQAAUe77op3Knpaj/xbt1MUD2+mFywerU4ukRnmv/xWc3o3y+gDQkI5YcKy1NcaYW1VbVlySXrfWLjfGTJK0xFo7S9ITkppLmmaMkaRCa+1F1tptxpiHVVuSJGnSvoEDAACg/vZ6a/SHT9z68+cFat0sXtN+cYIuH9JBgZ+/jcLrcSu2VRfFJCQ32nsAQEMJ6hoca+0cSXMOeu6BOl+fdZhjX5f0+tEGBAAAteavKdWEaTlaW7ZX1w/vqidG91dacnyjv2+Vx62EDn0b/X0AoCGExygEAACi2La9Xt01e4Ve/3ajerdupvk3naQzerdukve21srrcavFSVcfeWcACAEUHAAAQpS1VtNyivXLGXkq21ute/+vtx44J0NJca4my+ArL5V/7w5u8gkgbFBwAAAIQUU7KnTz9DzNXlGiEzq30NwJIzS0U4smz8EENQDhhoIDAEAI8futXv56g+79cKVq/H49OXqAfnVKD8W6YhzJs6/gJLCCAyBMUHAAAAgRK0vKdcPUHH25frvOzmitl68Yop7pzRzNVOVxS644xbXu5mgOAAgWBQcAAId5a/yaPD9ff/x0jZonuPT/fjpUPz+hc6OOfg46m2e14tv2knHxKwOA8MB3KwAAHPT1+m26YVqulnvK9dPjOumZiweqbUqC07H283rcDBgAEFYoOAAAOKC8ska//2iVnv9ynTq3SNQH1w/XqAHtnI51AOv3ybslX82HnO90FAAIGgUHAIAm9uGKEt00PVdFOyt168k99Mfz+yklMfR+JFeXbZStrmIFB0BYCb3vpgAARKgt5VX61cxlenvpZg1sn6KvfpGlEd3SnI71o/aPiKbgAAgjFBwAABqZtVb/WFKkO2ct1+4qnyad11f3nNFb8bHOjH4O1v9GRPd1OAkABI+CAwBAIyoo26Mbp+Xq0zWlOrl7ml4Zm6n+7VKcjhUUb4lbMYkpcrUIrWuDAOBwKDgAADSCGp9ff1m4Tvd/vEqxMTF68fLBunFEN8XEOD/6OVhVgQlqoTCuGgCCRcEBAKCBLd20U9lTc/Rd0U5dNLCdXrhssDq3THI6Vr15PW4l9RrhdAwAqBcKDgAADaSi2qc/zHXryc/XqnWzeL3z8xM0JrNDWK6A+KurVF26Xi1O/oXTUQCgXig4AAA0gAX5pZowLVf5pXt0/fCuemJ0f6Ulxzsd66hVb1krWasEJqgBCDMUHAAAjsH2vV7dNXulXvu2UL3Sk/XviSfp//q0djrWMatiRDSAMEXBAQDgKFhrNT23WLfOWKbSPV7dc0ZvPXhuhpLiXE5HaxBez2pJUny7Pg4nAYD6oeAAAFBPm3ZW6JbpeXp/eYmO79xCH2WfqOM6t3A6VoPyetxytWgnV3JkfS4AkY+CAwBAkPx+qynfbNA9H65Utc+vJy4coNtP7aFYV2jfsPNoVHncSmjH6WkAwg8FBwCAIKwqKdcN03K1aN02ndmntf56xRD1at3M6ViNxutxK2XoaKdjAEC9UXAAADgMb41fjy3I1yPz1qhZvEt/u3Koxg3rHJajn4Pl27NDvl1bGDAAICxRcAAA+BHfbNiu7Kk5Wu4p11VDO+qZSwapXUqC07EanbdkjSQpoUNfh5MAQP1RcAAAOEh5ZY3u+3iVnlu0Tp1SEzX7+uG6cEA7p2M1GS8jogGEMQoOAAB1zFlZopum52njjgrdMrK7Hr2gv1ISo+vHZZXHLZkYxbXp6XQUAKi36PqODQDAj9hSXqXb31+uf32/Sf3bNdeiW07WyB6tnI7lCK/HrbjW3RUTF/mn4wGIPBQcAEBUs9bqje+KdMf7y1VeVaOHzsnQvWf2VkJsZNyw82h4PW5OTwMQtig4AICota5sr258N0fz3KU6qVuaXh2bqQHtU5yO5ShrrbwlbrXM+InTUQDgqFBwAABRp8bn17OL1un+j1crxkgvXDZYE0/qppiYyB39HKyanR75K3ezggMgbFFwAABRJWfzTmVPzdGSjTt14YB2evGyweqSluR0rJDhLV4tiQlqAMIXBQcAEBUqqn16eJ5bjy9Yq/TkOL19zfEaO7RjRN+w82gwIhpAuKPgAAAi3mf5pZowLVdrSvfoumFd9ORFA9QqOd7pWCGpyuOWiUtQXKsuTkcBgKNCwQEARKwdFdW6a/YKvfqfQvVMT9a8G0forIw2TscKaV6PW/Ht+sjExDgdBQCOCgUHABBxrLV6L69Yt763TFt2V+mu03vpoXMzlBzPj70j8Za4ldBpkNMxAOCo8Z0eABBRNu2s0K3vLdPMZR4d1ylVH2YP1/GdWzodKyxYX428W9Yq5YTLnI4CAEeNggMAiAh+v9Ur/9mguz9YKW+NX4+N6q87T+upWBenWgWrunS95KtRAgMGAIQxCg4AIOyt3rJbN0zL0cKCbTqjd7qmjMlU79bNnI4VdqqYoAYgAlBwAABhy1vj1xOf5WvSJ2uUHO/Sa2Mzdd3wLox+PkqMiAYQCSg4AICw9J8N25U9NUfLPOUam9lRf7lkoNqnJjodK6x5PW7FNEuTq3m601EA4KhRcAAAYWV3VY3u+2iVnl20Th1TE/X+dcN00aD2TseKCF6PWwntM1gBAxDWKDgAgLDx0coSTZyep8LtFbp5ZHf9aVQ/pSbGOR0rYlR53GrW73SnYwDAMaHgAABC3tbdVbrj/eV687+b1L9dcy269WSd3KOV07Eiir9qj2q2beT6GwBhj4IDAAhZ1lq9+d9Nun3mMu2qqtGD52Tot2f2VkKsy+loEcdbki+JAQMAwh8FBwAQktZv26uJ7+Zq7uqtGtEtTa+OzdTA9ilOx4pY+yaocQ8cAOGOggMACCk+v9WzCwt038erFWOk5y4dpJtGdpcrhgvfG9P/7oHTx+EkAHBsKDgAgJCRu3mXsqfmaPHGHRrVv61eunyIuqQlOR0rKng9bsW26qyYBG6QCiC8UXAAAI6rrPbp4XluPb5grdKS4/T2Ncdr7NCOjCtuQl6Pm+tvAEQECg4AwFGfry3VhGm5cm/do2uHddGTowcovVm807GijtfjVurwsU7HAIBjRsEBADhiR0W17v5ghV75plA9WiVr3o0jdFZGG6djRaWa3WXy7dnGCg6AiEDBAQA0uRl5xbrlvTyVlFfprtN76aFzM5Qcz48kp3j3Dxig4AAIf/w0AQA0mc07K/XLGXl6L8+joR1TNXv8cJ3QpaXTsaIeBQdAJKHgAAAand9v9ep/CnX3BytUVePXY6P6647TeirOFeN0NChQcFyxim/d3ekoAHDMKDgAgEa1estuTZiWoy8KtumM3umaMiZTvVszijiUVHncim/TUyY2zukoAHDMKDgAgEZR7fPriQVrNWmeW0lxLr06NlPjh3dh9HMI8hav5vQ0ABGDggMAaHDfFm5X9tQc5RWXa0xmBz17ySC1T010OhYOwfr98pasUbNBZzsdBQAaBAUHANBgdlfV6P6PV+nZhevUITVRM68bposHtXc6Fg6jZnuRbHWlEge2m0QAACAASURBVFjBARAhKDgAgAbx8aotmvhurjZsr9DNI7vrT6P6KTWRazpCXdX+CWp9HU4CAA2DggMAOCalu6t0x6zl+ud3m9SvbXMtvGWkftIz3elYCBIjogFEGgoOAOCoWGv11n836fb3l2tnZbUeODtDvzurtxJiXU5HQz14PW6ZhGaKbdnB6SgA0CAoOACAelu/ba8mvpuruau3akS3NL0yZogGdUh1OhaOgtfjVkL7DKbbAYgYQd1hzRhznjFmtTEm3xhz7yG2n2qM+a8xpsYYc8VB23zGmKWBP7MaKjgAoOn5/FbPfFGggU98pi/Xb9Nzlw7SoltPptyEsSqPm9PTAESUI67gGGNckl6QdLakIkmLjTGzrLUr6uxWKOlaSb85xEtUWGuHNkBWAICDcjfvUvbUHC3euEMX9G+rly4frK5pyU7HwjGwNV5Vb12nFidd7XQUAGgwwZyiNlxSvrW2QJKMMW9LuljS/oJjrV0f2OZvhIwAAAdVVvv0yKdr9Nj8fKUlx+lf1xyvK4d25JSmCODdUiBZv+LbsYIDIHIEU3A6SdpY53GRpBPr8R6JxpglkmokTbbWzjx4B2PMBEkTJKlr1671eGkAQGP6Ym2ZbpiWI/fWPRqX1Vl/vmig0pvFOx0LDWTfBDXugQMgkjTFkIFu1tpNxpiekuYbY/KstWvr7mCtnSJpiiRlZWXZJsgEADiMnRXVuufDlfrr1xvUo1WyPpkwQmf3beN0LDSw/90Dp4/DSQCg4QRTcDZJ6lLncefAc0Gx1m4K/LfAGPOZpOMkrT3sQQAAx8zIK9Yt7+WppLxKvzm9lx46J0PNEhi6GYm8ntVypbSRq1ma01EAoMEE8xNrsaQ+xpgeqi02V0kK6mpEY0yapL3W2ipjTGtJJ0t6/GjDAgAaz+adlfrljDy9l+fR0I6pmj1+uE7o0tLpWGhEXiaoAYhARyw41toaY8ytkuZKckl63Vq73BgzSdISa+0sY8wwSTMkpUkabYz5g7V2oKT+kv4aGD4Qo9prcFb8yFsBABzg91u99m2h7pq9QlU1fk0e1V93ntZTca6g7iSAMOb1uNU88wKnYwBAgwrqnANr7RxJcw567oE6Xy9W7alrBx/3laTBx5gRANBI3Ft3a8K0XH2+tkyn90rXlDFD1KdNc6djoQn4KnapZqeHFRwAEYeTqgEgClX7/HpiwVpNmudWUpxLr47N1PjhXRj9HEW8njWSxIhoABGHggMAUWZx4Q5lT81RbvEujcnsoGcvGaT2qYlOx0IT8+6foEbBARBZKDgAECX2VNXo/o9X6y8LC9Q+JVEzrxumiwe1dzoWHOL1uCVjFN+2l9NRAKBBUXAAIArMXbVFE6fnav22Ct00spv+dEF/tUiKczoWHFRV4lZcejfFxLN6ByCyUHAAIIKV7q7SnbNW6I3vitSvbXN9cctIndIz3elYCAGMiAYQqSg4ABCBrLX61/eb9KuZy7Wzslr3n91HvzuzjxLjXE5HQwiw1srrcavFyb9wOgoANDgKDgBEmA3b9uqm6Xn6aNUWndi1pV4dm6lBHVKdjoUQ4tu1Rf6KXUpgBQdABKLgAECE8Pmtnl+0Tr//aJUk6dlLBunmk7vLFcPoZxyoiglqACIYBQcAIkBe8S5lT83Rt4U7dH6/tnrp8sHq1irZ6VgIUd7i1ZIoOAAiEwUHAMJYZbVPf/x0jSbPz1fLpDi99bPjddVxHblhJw7LW+KWiUtQXHpXp6MAQIOj4ABAmFpYUKYbpuZo9dY9+kVWZ/159AC1bp7gdCyEAa/Hrfi2vWViGDoBIPJQcAAgzOysqNY9H67UX7/eoO6tkjR3wok6p29bp2MhjFR53Ero0M/pGADQKCg4ABBGZuYV65b3lslTXqk7T+upSef2VbMEvpUjeNbvk7ckXynHXeR0FABoFPxUBIAwULyrUr+csUzTc4s1pEOqZl43TMO6tnQ6FsJQdekGyVfNgAEAEYuCAwAhzFqr1/5TqN/MXqHKGr8evaCffnN6L8W5YpyOhjDlDYyI5h44ACIVBQcAQtSarbs1YVquPltbptN6pWvKmCHKaNPc6VgIc9wDB0Cko+AAQIip9vn158/W6qFP3EqMjdErY4Zo/PCuiuGGnWgAXo9bMckt5Epp43QUAGgUFBwACCFLNu5Q9tQc5WzepcuHdNBzlw5Sh9REp2Mhgng9bsW3y+BeSQAiFgUHAELAnqoaPTB3tZ75okDtUhL03rVZunRwB6djIQJ5S9xK7vMTp2MAQKOh4ACAwz5ZvUU3vpur9dsqdONJ3TR5VH+1TIpzOhYikN9boeqyQsWfwvU3ACIXBQcAHFK2x6s7Zy3XP5YUKaNNM31+80id2ivd6ViIYN6SfMlaBgwAiGgUHABoYtZa/ev7Tbr9/eXavrda953VR78/q48S41xOR0OE2zciOr5DX4eTAEDjoeAAQBMq3L5XN03P05yVWzS8a0v9e2KmBndIdToWosT+e+C06+NwEgBoPBQcAGgCPr/VC1+u0+/mrJKV9MzFA3XrT3rIxehnNKEqj1uxLTsqJpH7KQGIXBQcAGhky4p3KXtqjv5TuEPn9m2jl68You6tkp2OhSjk9bi5/gZAxKPgAEAjqarx6Y+frtHk+flKTYjVP68+Tlcf34n7j8Ax3hK3Uk64zOkYANCoKDgA0AgWFZTphmm5WrVlt645oZOeumig2jRPcDoWophv9zb5ykuVwAoOgAhHwQGABrSzolq/nbNSL321Qd3SkvTxDSfq3H5tnY4FqKpkjSRxihqAiEfBAYAG8v4yj255L0/Fuyp1x6k9Nem8vmqewLdZhIb9I6IpOAAiHD95AeAYeXZV6pczlund3GIN7pCi967N0vCuaU7HAg7g9bilGJfi2/RwOgoANCoKDgAcJWutXv92o34ze4Uqqn364/n9dNcZvRTninE6GvADXo9b8W16yMTGOx0FABoVBQcAjkJ+6R5NmJajBfllOrVnK00Zk6m+bbm3CEIXI6IBRAsKDgDUQ7XPr6c+L9BDc1crITZGU8YM0fXDuyqGG3YihFm/X1Uet5L7n+F0FABodBQcAAjSdxt3KHtqjpZu3qXLBrfXc5cOVscWiU7HAo6oZsdmWe9eVnAARAUKDgAcwV5vjR6c69ZTn69Vu5QETR+XpcuGdHA6FhC0fRPUuAcOgGhAwQGAw5i3eqtufDdX67bt1YQRXfXYhQPUMinO6VhAvVQxIhpAFKHgAMAhlO3x6tezluv/LSlSRptm+vzmkTq1V7rTsYCj4vW4ZeKTFduyo9NRAKDRUXAAoA5rrd5Zulm3zVym7Xur9fuz+ui+s/ooMc7ldDTgqNVOUOsjE8MIcwCRj4IDAAGF2/fq5ul5+nDlFg3r0lKf3pipIR1TnY4FHDOvx63ErkOdjgEATYKCAyDq+fxWL321Xr+ds1J+Kz198UD98ic95GL0MyKAramWd2uBUoePdToKADQJCg6AqLbcU64bpubo6w3bdW7fNnr5iiHq3irZ6VhAg/GWrpP8PgYMAIgaFBwAUamqxqdHP83Xn+avUWpCrN64+jj97PhOMoZVG0QWLxPUAEQZCg6AqPPlum26YVqOVpbs1jUndNJTFw1Um+YJTscCGgX3wAEQbSg4AKLGrspq/fbDVXrxq/Xqlpakj244Uef1a+t0LKBReT1uuVJay9W8ldNRAKBJUHAARIVZyzy6+b08bd5VqdtP7aGHz+un5gl8C0Tkqyperfh2rN4AiB78dAcQ0Ty7KnXbzGWallOswR1SNH1clk7sluZ0LKDJeD1uNRt0jtMxAKDJUHAARCRrrf727Ub9evYKVVT79Mfz++muM3opzsWNDhE9/JW7VbNjM9ffAIgqFBwAEWdt6R5NmJar+fmlOqVnK70yJlN92zZ3OhbQ5KpK1khighqA6ELBARAxanx+PfV5gR6cu1rxsTF6+YrBuuHEborhhp2IUoyIBhCNKDgAIsJ/i3Yoe2qOvt+0S5cMaq/nLxukTi2SnI4FOGp/wWnX2+EkANB0KDgAwtpeb40emuvWU18UqE2zeE0fl6XLhnRwOhYQErwet+LSuyomnrIPIHpQcACErX+7t2rCu7kqKNurG0Z01eMXDlDLpDinYwEho8rj5vQ0AFGHggMg7Gzb69WvZ63Q3xdvVJ/WzbTgppN0eu/WTscCQoq1Vl6PWy1OutrpKADQpCg4AMKGtVZTl27WbTOXadveav3uzN667+wMJcW5nI4GhBxfean8e3ewggMg6lBwAISFjdsrdPN7efpgRYmyurTQJzeOUGbHFk7HAkLW/wYMUHAARBcKDoCQ5vdbvfTVet07Z6X8VvrzRQN02096KJYbdgKHta/gJHTo63ASAGhaFBwAIWu5p1w3TM3R1xu265yMNnr5iiHqkZ7sdCwgLFR5VkuuOMW17uZ0FABoUhQcACGnqsanP/07X4/+e41SEmL1j58O1TUndJYx3LATCJbX41Z8u94yMVyjBiC6UHAAhJSv1m1T9rQcrSzZrauP66SnLx6otikJTscCwo7X41YCAwYARCEKDoCQsKuyWr+bs0ovfrVeXVomaU72cJ3fv53TsYCwZP0+ebfkq3nmBU5HAYAmF9RVusaY84wxq40x+caYew+x/VRjzH+NMTXGmCsO2jbOGLMm8GdcQwUHEDlmL/dowOOf6cWv1uu2n/TQ8rtOp9wAx6C6bKNsdRUT1ABEpSOu4BhjXJJekHS2pCJJi40xs6y1K+rsVijpWkm/OejYVpIelJQlyUr6LnDs9oaJDyCclZRX6bYZyzQ1Z7MGtU/R9HFZOrFbmtOxgLC3f0Q0p6gBiELBnKI2XFK+tbZAkowxb0u6WNL+gmOtXR/Y5j/o2HMlzbPWbgtsnyfpPEn/OubkTWDr+w+reus6p2MAEcfKas3WPVq8cYcyfVbjOqVqcItUueYZbXY6HBABqkrWSBLX4ACISsEUnE6SNtZ5XCTpxCBf/1DHdjp4J2PMBEkTJKlr165BvnTjqyj4VpWFOU7HACJKjd9q295q+Wp8Oik2Rq2S4xW73aiCdV2gQTUbeLZcLTjVE0D0CYkhA9baKZKmSFJWVpZ1OM5+Xe+Y7XQEIGLU+Px6+osCPTh3tWKbxejxC/trwohuiolh9DMAAGg4wRScTZK61HncOfBcMDZJOv2gYz8L8lgAEeK/RTuUPTVH32/apYsHttMLlw9WpxZJTscCAAARKJiCs1hSH2NMD9UWlqskXR3k68+V9KgxZt9Vw+dI+m29UwIIS3u9NXporltPfVGgNs3i9e64E3TZ4A7csBMAADSaIxYca22NMeZW1ZYVl6TXrbXLjTGTJC2x1s4yxgyTNENSmqTRxpg/WGsHWmu3GWMeVm1JkqRJ+wYOAIhs89eUasK0HK0t26vsE7vq8Qv7Ky053ulYAAAgwhlrQ+aSF0m11+AsWbLE6RgAjtK2vV7dNXuFXv92o3q3bqYpY4bojN6tnY4FAAAijDHmO2tt1sHPh8SQAQDhz1qraTnF+uWMPJXtrda9/9dbD5yToaQ4l9PRAABAFKHgADhmRTsqdPP0PM1eUaITOrfQ3AkjNLRTC6djAQCAKETBAXDU/H6rl7/eoHs/XKkav19Pjh6gX53SQ7GuGKejAQCAKEXBAXBUVpaUK3tqjr5av11nZ7TWy1cMUc/0Zk7HAgAAUY6CA6BevDV+TZ6frz9+ukbNE1z6fz8dqp+f0JnRzwAAICRQcAAE7ev123TDtFwt95Trp8d10jMXD1TblASnYwEAAOxHwQFwROWVNfr9R6v0/Jfr1LlFoj64frhGDWjndCwAAIAfoOAAOKwPV5Topum5KtpZqVtP7qE/nt9PKYl86wAAAKGJ31IAHNKW8ir9auYyvb10swa2T9FXv8jSiG5pTscCAAA4LAoOgANYa/WPJUW6c9Zy7a7yadJ5fXXPGb0VH8voZwAAEPooOAD2Kyjboxun5erTNaU6uXuaXhmbqf7tUpyOBQAAEDQKDgDV+Pz6y8J1uv/jVYqNidGLlw/WjSO6KSaG0c8AACC8UHCAKLd0005lT83Rd0U7NXpAO714+WB1bpnkdCwAAICjQsEBolRFtU+TPnHric/WqnWzeE39xQm6YkgHbtgJAADCGgUHiEIL8ks1YVqu8kv3aPzwLnpi9AC1So53OhYAAMAxo+AAUWT7Xq/umr1Sr31bqF7pyfr3xJP0f31aOx0LAACgwVBwgChgrdX03GLdOmOZSvd4dc8ZvfXguRlKinM5HQ0AAKBBUXCACLdpZ4VumZ6n95eX6PjOLfRR9ok6rnMLp2MBAAA0CgoOEKH8fqsp32zQPR+uVLXPrycuHKDbT+2hWBc37AQAAJGLggNEoFUl5bphWq4WrdumM/u01l+vGKJerZs5HQsAAKDRUXCACOKt8euxBfl6ZN4aNYt36W9XDtW4YZ0Z/QwAAKIGBQeIEN9s2K7sqTla7inXVUM76plLBqldSoLTsQAAAJoUBQcIc+WVNbrv41V6btE6dUpN1Ozrh+vCAe2cjgUAAOAICg4QxuasLNFN0/O0cUeFbhnZXY9e0F8pifxvDQAAohe/CQFhaEt5lW5/f7n+9f0mDWjXXF/eerJO6t7K6VgAAACOo+AAYcRaqze+K9Id7y9XeVWNHjonQ/ee2VsJsdywEwAAQKLgAGFjXdle3fhujua5S3VStzS9OjZTA9qnOB0LAAAgpFBwgBBX4/Pr2UXrdP/HqxVjpBcuG6yJJ3VTTAyjnwEAAA5GwQFC2NJNO5U9NUffFe3UhQPa6cXLBqtLWpLTsQAAAEIWBQcIQRXVPk36xK0nPlur9OQ4vfPzEzQmswM37AQAADgCCg4QYhbkl2rCtFzll+7RdcO66MmLBqhVcrzTsQAAAMICBQcIEdv3enX3Byv16n8K1TM9WZ/eOEJnZrRxOhYAAEBYoeAADrPW6r28Yt363jJt2V2lu07vpYfOzVByPP97AgAA1Be/QQEO2rSzQre+t0wzl3l0XKdUfZg9XMd3bul0LAAAgLBFwQEc4PdbTflmg+75cKW8NX49Nqq/7jytp2JdMU5HAwAACGsUHKCJrSop1w3TcrVo3Tb9X+/W+uuYIerdupnTsQAAACICBQdoIt4avx5fkK+H561Rs3iXXr8yU9cO68LoZwAAgAZEwQGawH82bFf21Bwt85RrbGZH/eWSgWqfmuh0LAAAgIhDwQEa0e6qGt330So9u2idOqUmatb4YRo9sL3TsQAAACIWBQdoJB+tLNHE6XnauKNCN4/srkcv6KfUxDinYwEAAEQ0Cg7QwLburtId7y/Xm//dpP7tmmvhLSfr5B6tnI4FAAAQFSg4QAOx1uqf3xXpjveXa1dVjR48J0O/PbO3EmJdTkcDAACIGhQcoAGsK9urie/m6hP3Vo3olqZXx2ZqYPsUp2MBAABEHQoOcAx8fqtnFxbovo9XK8ZIz106SDeN7C5XDKOfAQAAnEDBAY5S7uZdyp6ao8Ubd2hU/7Z66fIh6pKW5HQsAACAqEbBAeqpstqnh+e59fiCtUpLjtPb1xyvsUM7csNOAACAEEDBAerh87WlmjAtV+6te3TtsC56cvQApTeLdzoWAAAAAig4QBB2VFTr7g9W6JVvCtWjVbLm3ThCZ2W0cToWAAAADkLBAY5gRl6xbnkvTyXlVbrr9F566NwMJcfzvw4AAEAo4rc04Eds3lmpX87I03t5Hg3tmKrZ44frhC4tnY4FAACAw6DgAAfx+61e/U+h7v5ghapq/Jo8qr/uPK2n4lwxTkcDAADAEVBwgDpWb9mtCdNy9EXBNp3RO11TxmSqd+tmTscCAABAkCg4gKRqn19PLFirSfPcSopz6bWxmbpueBdGPwMAAIQZCg6i3reF25U9NUd5xeUak9lBz14ySO1TE52OBQAAgKNAwUHU2lNVo/s/Xq2/LCxQh9REzbxumC4e1N7pWAAAADgGFBxEpbmrtmji9Fyt31ahm0d2159G9VNqYpzTsQAAAHCMKDiIKqW7q3TnrBV647si9WvbXAtvGamf9Ex3OhYAAAAaCAUHUcFaq7f+u0m3v79cOyur9cDZGfrdWb2VEOtyOhoAAAAaEAUHEW/Dtr2aOD1XH6/aqhHd0vTKmCEa1CHV6VgAAABoBBQcRCyf3+r5Rev0+49WyRjpuUsH6aaR3eWKYfQzAABApAqq4BhjzpP0F0kuSa9aaycftD1B0j8knSCpTNKV1tr1xpjuklZKWh3Y9Rtr7cSGiQ78uLziXcqemqNvC3fogv5t9dLlg9U1LdnpWAAAAGhkRyw4xhiXpBcknS2pSNJiY8wsa+2KOrtdL2m7tba3MeYqSY9JujKwba21dmgD5wYOqbLap0c+XaPH5ucrLTlOb/3seF11XEdu2AkAABAlglnBGS4p31pbIEnGmLclXSypbsG5WNJDga/flfS84TdKNLEv1pbphmk5cm/do3FZnfXniwYqvVm807EAAADQhIIpOJ0kbazzuEjSiT+2j7W2xhizU9K+2bs9jDHfS9ol6T5r7cKD38AYM0HSBEnq2rVrvT4AsLOiWvd8uFJ//XqDerRK1icTRujsvm2cjgUAAAAHNPaQgWJJXa21ZcaYEyTNNMYMtNbuqruTtXaKpCmSlJWVZRs5EyLIjLxi3fJenkrKq/Tr03rqD+f2VbMEZmcAAABEq2B+E9wkqUudx50Dzx1qnyJjTKykFpLKrLVWUpUkWWu/M8aslZQhacmxBkd027yzUr+ckaf38jzK7JiqWeOHK6tLS6djAQAAwGHBFJzFkvoYY3qotshcJenqg/aZJWmcpK8lXSFpvrXWGmPaSNpmrfUZY3pK6iOpoMHSI+r4/VavfVuou2avUFWNX5NH9dedp/VUnCvG6WgAAAAIAUcsOIFram6VNFe1Y6Jft9YuN8ZMkrTEWjtL0muS3jDG5EvaptoSJEmnSppkjKmW5Jc00Vq7rTE+CCKfe+tuTZiWq8/Xlun0XumaMmaI+rRp7nQsAAAAhBBTexZZ6MjKyrJLlnAGG/6n2ufXEwvWatI8t5LiXHpy9ACNH96F0c8AAABRzBjznbU26+DnuRobIW1x4Q5lT81RbvEujcnsoGcvGaT2qYlOxwIAAECIouAgJO2pqtH9H6/WXxYWqH1KomZeN0wXD2rvdCwAAACEOAoOQs7cVVs0cXqu1m+r0E0ju+lPF/RXi6Q4p2MBAAAgDFBwEDJKd1fpzlkr9MZ3RerXtrkW3jJSP+mZfuQDAQAAgAAKDhxnrdW/vt+kX81crh0V1br/7D763Zl9lBjncjoaAAAAwgwFB47asG2vbpqep49WbdGJXVvq1bGZGtQh1elYAAAACFMUHDjC57d6ftE6/f6jVZKkZy8ZpJtP7i5XDKOfAQAAcPQoOGhyecW7lD01R98W7tD5/drqpcsHq1urZKdjAfj/7d17dJx1ncfxzzfJZDKTzKRNm8v0fqEXCi23cBUQXERwXUAEZfFSV2hRQdRFdxE8K15Q3LMWduW2uHIQFtTScikCghQXOIICaltKy6WUQoEkvUJuM8lcfvvHPJlOMlNbIM2TTN6vczgzeWaenO/Q5zydT3+/7+8HAEAJIOBgyCSSaV35yMu66tENGhMK6I5PH6pzDpnAhp0AAAAYNAQcDIknNm7XoqWr9eLWLn2ueZJ+8g/zNL4m6HdZAAAAKDEEHOxT78ST+tf71+u/n3pN0+pCemjxkTp5ToPfZQEAAKBEEXCwz9zzXIsuvGutWjsSuuSDM/Tdj8xRdZBLDgAAAPsO3zYx6FraE/rK3Wu1fE2LDpoQ1b1fOFzNk8f4XRYAAABGAQIOBo1zTj//0+v6xn3rlEhl9KOPztUlJ8xUoLzM79IAAAAwShBwMChe2tqpxXeu0WOvbNcHZ47TTWcv0Oz6Gr/LAgAAwChDwMH7kkxn9B//94q++/BLqqoo08/OXqAvHDFFZWzYCQAAAB8QcPCePfP62zp/6WqtaWnXJxbE9NOPH6hYtMrvsgAAADCKEXDwrnX1pPRvD72oax7fqKZIle7+fLPOmB/zuywAAACAgIN35+EXt+iCZWu0aUdcXzx6qq76+/1VGwr4XRYAAAAgiYCDvbS9q1f/vOJ53frsG5pTX63HLzxGx80Y53dZAAAAQD8EHPxNzjn98q9v6mv3Pq+d3Ul9+6RZuvykWaoKlPtdGgAAAFCAgIPdem1Ht760/Dk9+MIWHTFljFZ+8SDNj0X9LgsAIKk3nVJrvEOt8Q61dLerJd7uPc8e29nb7XeJAEpMfVWNln9ood9l7BEBBwXSGafr/vCqLnvgBUnSNacfoIuOna5yln4GgH3KOaf2ZEIt3R15gcV7jLd74aVdLfEO7egpDDAm0/iqsGKhqOqCYXHXBjCYym1k3FUIOOhnbUu7zl+6Wn96/W2dMrdeN3xigabVhf0uCwBGtFQmrS2JTrX2Cy7FQ0winSo4P1heoVgooqZQVLOj9Tq+aYZioahi4aiaQhHFQtnHhlCNAmVMIQYwuhFwIElKJNP64cqX9aOVGzQmFND/nnuIzj10omyEJHUA8ENXskctRaeJ9QWW7Gtb4p1ycgXnj60MKRaOKhaK6piGaYqFI/0CS99rtZVV3I8BYC8RcKAnNm7XoqWr9eLWLn32sElacto8ja8J+l0WAPgi4zLa3tPtBZYOtfY9Ful16Uz1FJxfYWVqCkXUFI5ocnWtjhg/ud9IS1+IaQpFFSznr2EAGGzcWUexd+JJXXr/et341GuaOjak3y46Uh+Z2+B3WQCwTyRSyV1N+XkjLdnw0p5rzm+LdyjlMgXnRwLB3MjKoeMmKTYpkhtlyY26hCMaFwyrzMp8+IQAAImAM2rdu7ZVF971nFraE/r68TP0vVPmqCbI5QBgZHHO6e3eeDawdGenhOX3s+RPF9vZGy8432RqCNUoZJaH7QAAEfxJREFU5gWUBXWxvOlhu6aKNYUiqg4wsg0AIwHfaEeZ1vaEvnL3Wi1b06L5sYju+nyzjpgy1u+yAKCfZCatLfFOb+WwgYGl/yhMT5Gm/FB5IDeysv+YRn0otl+unyV/1KW+qloVNOUDQEkh4IwSzjnd/PRmfeO+dYon07ry1Ln65okzFShnGgWAodPhLYGcPyWsf4jJvrYt0V20KX9cMJybCnZ8dEZ2ZbFw1FthbFdTfiQQpCkfAEYpAs4o8PLWTl2wbI1+v2G7jp9Rp5vOPkhzGmr8LgtAiUhnMtrW01Vkv5bCaWJdqd6C8wNl5V4PS0QzInU6pmHqgFXEsg35jaEaVdKUDwDYA/6mKGHJdEZLHtuoKx56UcGKMt109gKdd8QUlbFhJ4C9EE8lcyMtu990skNbEp1KF2nKr62sygWVw72VxPrCSn5/S10wzGgLAGDQEHBK1J83v63zl67Wqrfadeb8Jv304/M1obbK77IA+Mw5px093QP6WdqLNue/05soOL/MTI1VXgN+OKJDx00s2GwyuwxyVKGKgA+fEAAw2hFwSkx3b0rfeeglLXnsFTVGglq+sFlnLoj5XRaAfaw3nVKb15RfbL+WvhDTGu9QMpMuOD9cEfD2aIlqwdiYTp4wp+imk+OD1Sovo3cPADB8EXBKyO9e3KoLlq3Rqzu6tfioKfrxx+ZpTIh/QQVGKuec2pOJAfu1DBhp8aaPbe/pLvo76quqcwFlbm3DgJXEdu3dEgkwwgsAKA0EnBKwvatXl6x4Xr949g3Nrq/WY18+RsfPHOd3WQB2I53JaEuic7f7teRPF4unkwXnB8srciMr+0XH67im6QXLH8dCUTWEahRgCWQAwChDwBnBnHP69aq3dPE9a7WzO6nLT5qlb580S1UBvtAAfuhK9vRb6nh3yyBvTXQq4wqXQB5bGcqFlGMapvXrZ8kuh5wNLmMqQzTlAwCwGwScEer1nd368vLndP/6LTp88hg9csFBWjAh6ndZQMnJuIy293TnpoL136+l/3SxjmRPwfnlVqYmb4+WSdW1ah4/WbEBIy2xcESNVRFV0ZQPAMD7RsAZYdIZpxue3KRvPbBeGSddffoB+sqx01XO0s/Au9KTThWMshRrzm+LdyhVZAnkmopgrgn/4LoJOjU8N7eXS9/Gk7FwVOOCYZUZTfkAAAwVAs4I8nxrhxYtXa2nXtupj8yp141nLdC0urDfZQHDhnNOb/fGd7tfS36I2dkbLzjfZNmmfG905cAxTQV7tvS9VhMI+vAJAQDAnhBwRoCeVFo/fGSDfvToy4oGK3TbuYfo04dOZA4+Ro1UJq22eGeRfpa+VcR2LYPck04VnB8sr8iOqISimhOt1wlNMwsa8mPhiOqraMoHAGCkI+AMc394dYcW3bla69s69ZnDJmrJaQeovoZ/OUZp6Ez2FNmvpXC62LZEt5wKm/LrguFcQPlAw/SiIy2xUFS1lVX8gwAAAKMEAWeYak8k9a37X9D1T27S1LEhPbjoSJ0yt8HvsoA9yriMtia6drtfS0u8I/e8K9VbcH5FX1N+OKKpNWN1VP1Ub/WwSG66WJP3PFjOLQwAAPTHt4Nh6L7nW/Xl5c/pzfaEvnb8dH3/lLmqCfJHBX/FU0m15o20FO7dkn1si3cqXaQpPxII5qaCHTZ+khdYBk4Ti6ouGKIpHwAAvGd8ax5G2jp6dPHda7V09VuaH4to2cJmHTl1rN9loYQ557SzNz5gmljeMsjd2RGXlni73ulNFJxfZqaGqppcQDmobkI2qOTt2dK3RHI1TfkAAGAIEHCGAeecbnlmsy5ZsU5dvWn94NQ5+uYJ+6mygn/FxnvTm055TfntRfdr6XveGu9QbyZdcH6oPJAbWTlgTKNOmjArF2LyR13qq6pVQVM+AAAYRgg4PntlW5cW37lGj27YpuNm1OmmsxZobmPE77IwDDnn1JHsGdCAv2uUJT/EbOvpKvo7xgXDualgc2ob+u3Z0hdcYqGoIoEgTfkAAGBEIuD4JJXOaMljG/Wdh15UZUWZbjxrvhYdOVVlbNg56qQzGW1NdGangg3oZ8mfLtbS3a54OllwfqCsPBdQZkbH6QON0/qNtPQ15zeGalRJUz4AAChxfNvxwV/eeFvnL12tv77ZrjMObNK1Zx6oibUhv8vCIOtO9fZb6njXfi0d/XpetiQ6lXGFSyDXVlblgsqR9VO8KWL9m/KbQhHVBcOMtgAAAHgIOEOouzelKx56SUse36j66kotX9isMxfE/C4L74JzTjt6ur0+lo7iG096q4y1J4s35TdWRRQLRzQxXKvDxk0qCCzZZZCjClUEfPiEAAAAIxsBZ4isfGmrFi9bo43bu7XoqCn694/N05gQX2CHi76m/N2Flb7nrfEOJYs05VdXVGYDSjiiBWNjOnnCnKKbTo4PVqu8jMUjAAAA9hUCzj62o7tXl6xYp1ue2axZ46v1+y8drRP2G+93WaOCc07tyYQ30vK3VxPb3tNd9HfUV1XnQsq8MY0DRlkiuVATCVQN8acDAABAMQScfcQ5p6Wr3tLF96zVju6kLvu7/fTtD89WKMCSuu9XKpPW1kTXgKBSPMQUa8oPllfkwsnsaL2Ob5pRMNISC0XVEKpRgCWQAQAARhQCzj6weWdcF971nO5b16bmybV6+IKjdNCEWr/LGva6kj27Npj0Qkux1cS27qYpf2xlKBdQjmmY1q+fJb85f0xliKZ8AACAEkXAGUSZjNMNT27SpQ+sV8ZJS06bp4uPm6HyUbz0c8ZltL2ne1dAyQsqA1cT60j2FJxfYWVq9ILKpOpaHT5+cr/9WvpCTGNVRFU05QMAAIx6BJxBsq61Q4vuXK0nN+3UybPrdeNZCzR9XNjvsvaZRCqptkTHHlcSa4t3KOUyBedHAsHcVLBD6ibm7deSN00sHNW4YFhlRlM+AAAA9g4B533qSaV11coNunLly4oGK3TbuYfo04dOHJFToJxzers3njeysmuq2MDNJ3f2xgvON5kaQjW5oDJ/bGzA8se7Qkx1IOjDJwQAAECpI+C8D0++ukOL7lytdW2dOveQibrmjANUXzP8vrinMuncEsgDVxHr23yy77WedKrg/KryCsVCUcXCUc2tbdCJTTOzfS0DlkGur6pWBU35AAAA8BEB5z1oTyR12QMv6PonN2nymJAeOP8Inbp/45DX0ZFMeAGlcL+W/BCzLdEtp8Km/HHBcC6cHBedngsx2VGXSC7ERANVI3JECgAAAKMPAedd+s26Nn1p2Rq92Z7QxcdO1w9Onaua4OD9b8y4jLYmuors15K3DLIXZLpSvQXnB8rK1eRNA5tWU6ejG6YWmSaWXQI5WM4fPwAAAErLXn3DNbNTJP2npHJJ/+Ocu2rA60FJt0o6TNJ2SZ9yzm3yXvuWpPMkpSVd7Jx7aNCqH0JtHT366j1r9etVb+nApoiWLWzWkVPH7vX58VRyt/0s+SFmS6JT6SJN+dFAVW5KWPP4Sf32a8nfdHJsMERTPgAAAEatPQYcMyuXdJ2kD0t6Q9IzZrbCObcu723nSdrpnNvPzM6R9GNJnzKzeZLOkXSApAmSHjGz2c659GB/kH3FOadbntmsS1asU1dvWt8/ZY7+5cT9VFlRJuecdvR073GzyZZ4u97pTRT87jIzNVZFcksdH1w3sd9+LbFQVE1eeAlXVPrw6QEAAICRZW9GcI6QtME5t1GSzOxXkk6XlB9wTpd0hfd8maRrLdu0cbqkXznneiS9amYbvN/31OCUv2/95pVXdPnKVVqzZZumTC7XyTOr9Yz+oOMefDAXYpKZwqwWrgjk+lkOHNukD0+YnbfpZCT32vhgtcrLGG0BAAAABsveBJyJkjbn/fyGpCN39x7nXMrM3pE0zjv+xwHnTnzP1Q6hTMbprJW3qyfQLk2UXpcU31qdbbwPRbT/mMYBvS27mvIjgSq/ywcAAABGpWHRZW5miyUtlqQpU6b4XE1WWZnpuqPPVCxSpYPq69UQqlGAJZABAACAYW1vAs6bkibn/TzJO1bsPW+YWYWkWmUXG9ibc+Wcu0nSTZLU3NxcuJ6xT86bf6DfJQAAAAB4F/amAeQZSbPMbLqZVSq7aMCKAe9ZIWmh9/wsSY8655x3/BwzC5rZdEmzJD09OKUDAAAAQH97HMHxemoukvSQsstE3+yce97MvifpWefcCkk/l3Sbt4jADmVDkLz3LVV2QYKUpAtH0gpqAAAAAEYWyw60DB/Nzc3u2Wef9bsMAAAAAMOYmf3ZOdc88DhrFAMAAAAoGQQcAAAAACWDgAMAAACgZBBwAAAAAJQMAg4AAACAkkHAAQAAAFAyCDgAAAAASgYBBwAAAEDJIOAAAAAAKBkEHAAAAAAlg4ADAAAAoGQQcAAAAACUDHPO+V1DP2a2VdJrfteRZ7ykbX4XgZLHdYahwrWGocB1hqHCtTa6TXXO1Q88OOwCznBjZs8655r9rgOljesMQ4VrDUOB6wxDhWsNxTBFDQAAAEDJIOAAAAAAKBkEnD27ye8CMCpwnWGocK1hKHCdYahwraEAPTgAAAAASgYjOAAAAABKBgEHAAAAQMkg4OyGmZ1iZi+a2QYzu9TvelC6zGyTmT1nZqvM7Fm/60HpMLObzWyLma3NO1ZnZr8zs5e9x7F+1oiRbzfX2RVm9qZ3X1tlZh/1s0aMfGY22cx+b2brzOx5M/uqd5x7GgoQcIows3JJ10k6VdI8Sf9oZvP8rQol7kTn3MGs5Y9BdoukUwYcu1TSSufcLEkrvZ+B9+MWFV5nknS1d1872Dn3wBDXhNKTknSJc26epKMkXeh9N+OehgIEnOKOkLTBObfROdcr6VeSTve5JgB4V5xzj0vaMeDw6ZJ+4T3/haQzhrQolJzdXGfAoHLOtTjn/uI975C0XtJEcU9DEQSc4iZK2pz38xveMWBfcJIeNrM/m9liv4tByWt0zrV4z1slNfpZDEraRWa2xpvCxrQhDBozmybpEEl/Evc0FEHAAfx3rHPuUGWnRF5oZsf7XRBGB5fdJ4C9ArAv3CBppqSDJbVI+om/5aBUmFmNpOWSvuaca89/jXsa+hBwintT0uS8nyd5x4BB55x703vcIuluZadIAvtKm5nFJMl73OJzPShBzrk251zaOZeR9DNxX8MgMLOAsuHmdufcXd5h7mkoQMAp7hlJs8xsuplVSjpH0gqfa0IJMrNqM4v0PZd0sqS1f/ss4H1ZIWmh93yhpHt9rAUlqu8Lp+fj4r6G98nMTNLPJa13zi3Je4l7GgpYdjQPA3lLWl4jqVzSzc65K30uCSXIzGYoO2ojSRWS7uBaw2Axs19KOkHSeEltkr4j6R5JSyVNkfSapE8652gQx3u2m+vsBGWnpzlJmyRdkNcnAbxrZnaspCckPScp4x2+TNk+HO5p6IeAAwAAAKBkMEUNAAAAQMkg4AAAAAAoGQQcAAAAACWDgAMAAACgZBBwAAAAAJSMCr8LAACMTmaWVnbJ14CklKRbJV3tbQ4JAMB7QsABAPgl7pw7WJLMrEHSHZKiyu6jAgDAe8IUNQCA75xzWyQtlnSRZU0zsyfM7C/ef8dIkpndamZn9J1nZreb2elmdoCZPW1mq8xsjZnN8uuzAAD8xUafAABfmFmnc65mwLG3Jc2R1CEp45xLeGHll865ZjP7oKSvO+fOMLNaSaskzZJ0taQ/OuduN7NKSeXOufjQfiIAwHDAFDUAwHAUkHStmR0sKS1ptiQ55x4zs+vNrF7SJyQtd86lzOwpSZeb2SRJdznnXvatcgCAr5iiBgAYFsxshrJhZoukr0tqk3SQpGZJlXlvvVXSZyT9k6SbJck5d4ek0yTFJT1gZh8ausoBAMMJIzgAAN95IzI3SrrWOee86WdvOOcyZrZQUnne22+R9LSkVufcOu/8GZI2Ouf+y8ymSFog6dEh/RAAgGGBgAMA8EvIzFZp1zLRt0la4r12vaTlZvY5Sb+V1NV3knOuzczWS7on73d9UtJnzSwpqVXSD4egfgDAMMQiAwCAEcXMwsrun3Ooc+4dv+sBAAwv9OAAAEYMMztJ0npJPyXcAACKYQQHAAAAQMlgBAcAAABAySDgAAAAACgZBBwAAAAAJYOAAwAAAKBkEHAAAAAAlIz/BzJRNGffpoCBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[I 2020-03-08 22:08:38,416] Finished trial#0 resulted in value: 76476.9341496019. Current best value is 76476.9341496019 with parameters: {'isolation_rate': 0.27691874170288483, 'confirmation_rate': 0.02950017925660996}.\n",
      "[I 2020-03-08 22:12:30,333] Finished trial#1 resulted in value: 73281.080258802. Current best value is 73281.080258802 with parameters: {'isolation_rate': 0.26628882633274087, 'confirmation_rate': 0.0294126049399196}.\n",
      "[I 2020-03-08 22:15:50,246] Finished trial#2 resulted in value: 85265.64405348072. Current best value is 73281.080258802 with parameters: {'isolation_rate': 0.26628882633274087, 'confirmation_rate': 0.0294126049399196}.\n",
      "[I 2020-03-08 22:19:14,111] Finished trial#3 resulted in value: 79975.0195749794. Current best value is 73281.080258802 with parameters: {'isolation_rate': 0.26628882633274087, 'confirmation_rate': 0.0294126049399196}.\n",
      "[I 2020-03-08 22:22:31,003] Finished trial#4 resulted in value: 70869.31842937776. Current best value is 70869.31842937776 with parameters: {'isolation_rate': 0.25531803279599863, 'confirmation_rate': 0.029765754948419466}.\n",
      "[I 2020-03-08 22:25:58,054] Finished trial#5 resulted in value: 81624.82746831629. Current best value is 70869.31842937776 with parameters: {'isolation_rate': 0.25531803279599863, 'confirmation_rate': 0.029765754948419466}.\n",
      "[I 2020-03-08 22:29:34,529] Finished trial#6 resulted in value: 100639.01337222212. Current best value is 70869.31842937776 with parameters: {'isolation_rate': 0.25531803279599863, 'confirmation_rate': 0.029765754948419466}.\n",
      "[I 2020-03-08 22:32:55,571] Finished trial#7 resulted in value: 89237.61092050422. Current best value is 70869.31842937776 with parameters: {'isolation_rate': 0.25531803279599863, 'confirmation_rate': 0.029765754948419466}.\n",
      "[I 2020-03-08 22:36:34,225] Finished trial#8 resulted in value: 97019.26384150493. Current best value is 70869.31842937776 with parameters: {'isolation_rate': 0.25531803279599863, 'confirmation_rate': 0.029765754948419466}.\n",
      "[I 2020-03-08 22:40:02,336] Finished trial#9 resulted in value: 68886.26309723839. Current best value is 68886.26309723839 with parameters: {'isolation_rate': 0.25830960412786363, 'confirmation_rate': 0.029273275726341082}.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best_value =  68886.26309723839\n",
      "best_params =  {'isolation_rate': 0.25830960412786363, 'confirmation_rate': 0.029273275726341082}\n"
     ]
    }
   ],
   "source": [
    "#ignore warning messages\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "#import packages\n",
    "import optuna\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ddeint import ddeint\n",
    "import scipy.integrate as integrate\n",
    "from sklearn.metrics import mean_squared_error\n",
    "import matplotlib.pylab as plt\n",
    "import matplotlib.dates as mdates\n",
    "from matplotlib.pylab import rcParams\n",
    "plt.style.use('seaborn-colorblind')\n",
    "%matplotlib inline\n",
    "%load_ext autotime\n",
    "\n",
    "observation = 'COVID-19.csv'\n",
    "observed_J = pd.read_csv(observation)\n",
    "\n",
    "#parameters\n",
    "t_max = len(observed_J) - 1\n",
    "growth_rate = 0.31\n",
    "incidence_rate = 0.2\n",
    "confirmation_rate = 0.03\n",
    "isolation_rate = 0.25\n",
    "J0 = 21 #start day = 2020-02-14\n",
    "I0 = J0 * 10\n",
    "G0 = 0\n",
    "l1 = 0.1\n",
    "l2 = isolation_rate\n",
    "\n",
    "#set timespan\n",
    "tspan = np.linspace(0, t_max, t_max+1)\n",
    "\n",
    "#dynamic parameters\n",
    "def states_before_init(t):\n",
    "    factor = np.exp(-t**2)\n",
    "    return [I0 * factor, J0 * factor, G0 * factor]\n",
    "\n",
    "def f2(t):\n",
    "    return np.minimum(incidence_rate/14 * t, incidence_rate)\n",
    "\n",
    "def f4(t):\n",
    "    return np.minimum(confirmation_rate/14 * t, confirmation_rate)\n",
    "\n",
    "def lt(t):\n",
    "    return np.where(t < 15, l1, l2)\n",
    "\n",
    "f2map = map(lambda x: f2(x), tspan)\n",
    "f4map = map(lambda x: f4(x), tspan)\n",
    "ltmap = map(lambda x: lt(x), tspan)\n",
    "\n",
    "#plot dynamic parameters\n",
    "rcParams['figure.figsize'] = 14, 7\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(tspan, list(f2map))\n",
    "ax.plot(tspan, list(f4map))\n",
    "ax.plot(tspan, list(ltmap))\n",
    "ax.legend(['f2(t)', 'f4(t)', 'lt(t)'])\n",
    "ax.set_title('Dynamic Parameters')\n",
    "ax.set_xlabel('Days')\n",
    "plt.show()\n",
    "\n",
    "#set differential equation: Fudan CCDC model\n",
    "def Fudan_CCDC(v, t):\n",
    "    #I0 at time t must be positive\n",
    "    I0 = np.maximum(v(t)[0] - v(t)[1] - v(t)[2], 0)\n",
    "    #I0 at time s must be positive\n",
    "    Integral2, _ = integrate.quad(lambda s: f2(t - s) * np.maximum(v(s)[0] - v(s)[1] - v(s)[2], 0), 0, t)\n",
    "    Integral4, _ = integrate.quad(lambda s: f4(t - s) * np.maximum(v(s)[0] - v(s)[1] - v(s)[2], 0), 0, t)\n",
    "\n",
    "    dI = growth_rate * I0\n",
    "    dJ = growth_rate * Integral4\n",
    "    dG = lt(t) * (Integral2 - Integral4)\n",
    "    \n",
    "    return np.array([dI, dJ, dG])\n",
    "\n",
    "#fitting Fudan CCDC model to the Data\n",
    "def objective_FudanCCDC(trial):\n",
    "    global isolation_rate, confirmation_rate\n",
    "    isolation_rate = trial.suggest_uniform('isolation_rate', 0.1, 0.5)\n",
    "    confirmation_rate = trial.suggest_loguniform('confirmation_rate', 0.01, 0.1)\n",
    "    solution = ddeint(Fudan_CCDC, states_before_init, tspan)\n",
    "    return mean_squared_error(solution[:,1], observed_J['A'])\n",
    "\n",
    "study = optuna.create_study()\n",
    "study.optimize(objective_FudanCCDC, n_trials=10)\n",
    "print(\"best_value = \", study.best_value)\n",
    "print(\"best_params = \", study.best_params)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Days</th>\n",
       "      <th>Cumulative Infected</th>\n",
       "      <th>Observed Infected</th>\n",
       "      <th>Isolated Infected</th>\n",
       "      <th>Potential Infected</th>\n",
       "      <th>PCR Test Posives A</th>\n",
       "      <th>PCR Test Posives B</th>\n",
       "      <th>PCR Test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2020-02-14</td>\n",
       "      <td>210</td>\n",
       "      <td>21</td>\n",
       "      <td>0e+00</td>\n",
       "      <td>189</td>\n",
       "      <td>21</td>\n",
       "      <td>21</td>\n",
       "      <td>14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2020-02-15</td>\n",
       "      <td>269</td>\n",
       "      <td>21</td>\n",
       "      <td>4e-02</td>\n",
       "      <td>248</td>\n",
       "      <td>21</td>\n",
       "      <td>21</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2020-02-16</td>\n",
       "      <td>327</td>\n",
       "      <td>21</td>\n",
       "      <td>3e-01</td>\n",
       "      <td>306</td>\n",
       "      <td>21</td>\n",
       "      <td>21</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2020-02-17</td>\n",
       "      <td>386</td>\n",
       "      <td>22</td>\n",
       "      <td>1e+00</td>\n",
       "      <td>363</td>\n",
       "      <td>46</td>\n",
       "      <td>38</td>\n",
       "      <td>273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2020-02-18</td>\n",
       "      <td>444</td>\n",
       "      <td>22</td>\n",
       "      <td>2e+00</td>\n",
       "      <td>420</td>\n",
       "      <td>53</td>\n",
       "      <td>44</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2020-02-19</td>\n",
       "      <td>503</td>\n",
       "      <td>24</td>\n",
       "      <td>5e+00</td>\n",
       "      <td>475</td>\n",
       "      <td>60</td>\n",
       "      <td>50</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2020-02-20</td>\n",
       "      <td>601</td>\n",
       "      <td>26</td>\n",
       "      <td>9e+00</td>\n",
       "      <td>565</td>\n",
       "      <td>70</td>\n",
       "      <td>60</td>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2020-02-21</td>\n",
       "      <td>769</td>\n",
       "      <td>30</td>\n",
       "      <td>2e+01</td>\n",
       "      <td>721</td>\n",
       "      <td>79</td>\n",
       "      <td>69</td>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2020-02-22</td>\n",
       "      <td>976</td>\n",
       "      <td>36</td>\n",
       "      <td>3e+01</td>\n",
       "      <td>911</td>\n",
       "      <td>90</td>\n",
       "      <td>80</td>\n",
       "      <td>85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2020-02-23</td>\n",
       "      <td>1252</td>\n",
       "      <td>45</td>\n",
       "      <td>5e+01</td>\n",
       "      <td>1162</td>\n",
       "      <td>114</td>\n",
       "      <td>102</td>\n",
       "      <td>96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>2020-02-24</td>\n",
       "      <td>1607</td>\n",
       "      <td>57</td>\n",
       "      <td>7e+01</td>\n",
       "      <td>1483</td>\n",
       "      <td>126</td>\n",
       "      <td>114</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>2020-02-25</td>\n",
       "      <td>1967</td>\n",
       "      <td>72</td>\n",
       "      <td>1e+02</td>\n",
       "      <td>1799</td>\n",
       "      <td>140</td>\n",
       "      <td>128</td>\n",
       "      <td>104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>2020-02-26</td>\n",
       "      <td>2327</td>\n",
       "      <td>92</td>\n",
       "      <td>1e+02</td>\n",
       "      <td>2102</td>\n",
       "      <td>149</td>\n",
       "      <td>136</td>\n",
       "      <td>44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>2020-02-27</td>\n",
       "      <td>2846</td>\n",
       "      <td>118</td>\n",
       "      <td>2e+02</td>\n",
       "      <td>2545</td>\n",
       "      <td>171</td>\n",
       "      <td>156</td>\n",
       "      <td>168</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>2020-02-28</td>\n",
       "      <td>3606</td>\n",
       "      <td>153</td>\n",
       "      <td>2e+02</td>\n",
       "      <td>3206</td>\n",
       "      <td>195</td>\n",
       "      <td>180</td>\n",
       "      <td>151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>2020-02-29</td>\n",
       "      <td>4539</td>\n",
       "      <td>197</td>\n",
       "      <td>3e+02</td>\n",
       "      <td>4012</td>\n",
       "      <td>215</td>\n",
       "      <td>197</td>\n",
       "      <td>130</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>2020-03-01</td>\n",
       "      <td>5782</td>\n",
       "      <td>252</td>\n",
       "      <td>6e+02</td>\n",
       "      <td>4936</td>\n",
       "      <td>224</td>\n",
       "      <td>206</td>\n",
       "      <td>178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>2020-03-02</td>\n",
       "      <td>7278</td>\n",
       "      <td>323</td>\n",
       "      <td>9e+02</td>\n",
       "      <td>6030</td>\n",
       "      <td>239</td>\n",
       "      <td>221</td>\n",
       "      <td>96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>2020-03-03</td>\n",
       "      <td>9071</td>\n",
       "      <td>411</td>\n",
       "      <td>1e+03</td>\n",
       "      <td>7320</td>\n",
       "      <td>253</td>\n",
       "      <td>230</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>2020-03-04</td>\n",
       "      <td>11094</td>\n",
       "      <td>521</td>\n",
       "      <td>2e+03</td>\n",
       "      <td>8714</td>\n",
       "      <td>269</td>\n",
       "      <td>246</td>\n",
       "      <td>3835</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>2020-03-05</td>\n",
       "      <td>13635</td>\n",
       "      <td>659</td>\n",
       "      <td>3e+03</td>\n",
       "      <td>10473</td>\n",
       "      <td>302</td>\n",
       "      <td>276</td>\n",
       "      <td>258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>2020-03-06</td>\n",
       "      <td>16534</td>\n",
       "      <td>828</td>\n",
       "      <td>3e+03</td>\n",
       "      <td>12404</td>\n",
       "      <td>333</td>\n",
       "      <td>302</td>\n",
       "      <td>699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>2020-03-07</td>\n",
       "      <td>20154</td>\n",
       "      <td>1038</td>\n",
       "      <td>4e+03</td>\n",
       "      <td>14828</td>\n",
       "      <td>392</td>\n",
       "      <td>355</td>\n",
       "      <td>553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>2020-03-08</td>\n",
       "      <td>24701</td>\n",
       "      <td>1296</td>\n",
       "      <td>6e+03</td>\n",
       "      <td>17904</td>\n",
       "      <td>439</td>\n",
       "      <td>397</td>\n",
       "      <td>147</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>2020-03-09</td>\n",
       "      <td>29890</td>\n",
       "      <td>1612</td>\n",
       "      <td>7e+03</td>\n",
       "      <td>21292</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Days  Cumulative Infected  Observed Infected  Isolated Infected  \\\n",
       "0  2020-02-14                  210                 21              0e+00   \n",
       "1  2020-02-15                  269                 21              4e-02   \n",
       "2  2020-02-16                  327                 21              3e-01   \n",
       "3  2020-02-17                  386                 22              1e+00   \n",
       "4  2020-02-18                  444                 22              2e+00   \n",
       "5  2020-02-19                  503                 24              5e+00   \n",
       "6  2020-02-20                  601                 26              9e+00   \n",
       "7  2020-02-21                  769                 30              2e+01   \n",
       "8  2020-02-22                  976                 36              3e+01   \n",
       "9  2020-02-23                 1252                 45              5e+01   \n",
       "10 2020-02-24                 1607                 57              7e+01   \n",
       "11 2020-02-25                 1967                 72              1e+02   \n",
       "12 2020-02-26                 2327                 92              1e+02   \n",
       "13 2020-02-27                 2846                118              2e+02   \n",
       "14 2020-02-28                 3606                153              2e+02   \n",
       "15 2020-02-29                 4539                197              3e+02   \n",
       "16 2020-03-01                 5782                252              6e+02   \n",
       "17 2020-03-02                 7278                323              9e+02   \n",
       "18 2020-03-03                 9071                411              1e+03   \n",
       "19 2020-03-04                11094                521              2e+03   \n",
       "20 2020-03-05                13635                659              3e+03   \n",
       "21 2020-03-06                16534                828              3e+03   \n",
       "22 2020-03-07                20154               1038              4e+03   \n",
       "23 2020-03-08                24701               1296              6e+03   \n",
       "24 2020-03-09                29890               1612              7e+03   \n",
       "\n",
       "    Potential Infected  PCR Test Posives A  PCR Test Posives B  PCR Test  \n",
       "0                  189                  21                  21        14  \n",
       "1                  248                  21                  21         0  \n",
       "2                  306                  21                  21         0  \n",
       "3                  363                  46                  38       273  \n",
       "4                  420                  53                  44        36  \n",
       "5                  475                  60                  50        71  \n",
       "6                  565                  70                  60        90  \n",
       "7                  721                  79                  69        90  \n",
       "8                  911                  90                  80        85  \n",
       "9                 1162                 114                 102        96  \n",
       "10                1483                 126                 114        39  \n",
       "11                1799                 140                 128       104  \n",
       "12                2102                 149                 136        44  \n",
       "13                2545                 171                 156       168  \n",
       "14                3206                 195                 180       151  \n",
       "15                4012                 215                 197       130  \n",
       "16                4936                 224                 206       178  \n",
       "17                6030                 239                 221        96  \n",
       "18                7320                 253                 230        71  \n",
       "19                8714                 269                 246      3835  \n",
       "20               10473                 302                 276       258  \n",
       "21               12404                 333                 302       699  \n",
       "22               14828                 392                 355       553  \n",
       "23               17904                 439                 397       147  \n",
       "24               21292                 NaN                 NaN       NaN  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAG5CAYAAACjlISHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3iUZfr28e9FKKGEjtKLiAiBJJCAIKIgCthQkGIBxbW7+rOtBbu7urqvrr0tKCJiAVGxgbJWEKQkgHQFJECogdACIaTc7x8zZAOmDJDJk5mcn+PI4cxTzyE5JBf3c1+3OecQEREREREpTyp4HUBERERERKS0qRASEREREZFyR4WQiIiIiIiUOyqERERERESk3FEhJCIiIiIi5Y4KIRERERERKXdUCImISFgzszfM7GGvc4iISNmiQkhEpBwxsyvMLNHM0s1ss5lNM7Mz8u1vb2afm9luM9trZj+Y2en+fZFmtsvMzi7gus+b2WT/62QzO8f/eqSZ5fjvl25ma83sbTM7pZicQ81stpntN7MfC9h/kZkt9V9ztpm1L+xazrmbnHP/CPgP6fD7jDOzJ47lXBERKdtUCImIlBNmdhfwAvBP4ESgOfAacLF/f2tgFrAEaAU0Bj4FpptZd+fcAWAicNUR140ALgfeKeTWvzjnagC1gHOADCDJzDoUETfNn/XpAj5HG+A94CagNvAF8LmZVSzmj0BERCSPCiERkXLAzGoBfwf+6pz7xDm3zzmX5Zz7wjl3j/+wx/AVLQ8659Kcc3udcy8B7wL/8h/zDnCpmVXLd/l++P4+mVZUBudcjnNujXPuFuAn//0KO/Zb59wkYFMBu/sBM51zPzvnsv3ZmgBnFfLZ80Z1zKyXmaWY2d1mts0/KnZNUbmPuNZHZrbFP2I2w8yij7jPG2b2X/9o2k9m1iLf/hfNbIOZ7TGzJDPrmW/fY2Y2yczG+89dZmYJgeYSEZGjp0JIRKR86A5E4hvhKcy5wEcFbJ8E9DCzqs652cBmYFC+/SOA9/1FSaA+AXoWe1Th7IjXBhQ1wpRfQ3yjU02Aa4FXzaxOgOdOA9oAJwAL8I1M5Xcl8A+gPrDoiP3zgTigLvA+8JGZRebbPwD4EN8o1+fAKwFmEhGRYxCShZCZjfX/S97SAI5t7n/GfaGZLTaz80sjo4hIGVMP2F5MsVIfX5FzpM34/r6o638/Hv/jcWZWE9+jdYU9FleYTfmud7S+Bc7yj+5UBh4AKgPVij4tTxbwd/+I2FQgHWgbyInOubH+kbJMfCNasf7RtkO+cs7N8O9/EOhuZs38505wzu1wzmU75/4NVDnivj8756Y653LwjcLFBvh5RETkGIRkIQSMA/oHeOxDwCTnXCfgMnzPw4uIlDc7gPrFzKPZDjQqYHsjIBfY6X//LtDbzBoDg4E1zrmFR5mnCb55QIe6uh1qpvBAcSc651YCV+MbMdmMr4BbDqQEeO8dRxSE+4EaxZ1kZhFm9rSZrTGzPUCyf1f9fIdtyJczHd9nbOw//29mtsL/WN0ufKNS+c/dckSmSM17EhEJnpAshJxzM/D/BXqImbU2s6/9z13PNLNTDx0O1PS/rkXBz5uLiIS7X4BM4JIijvkWGFLA9qH45g7tB3DOrQNmAsPxPRZ3tKNBAAP91zjU1a2G/+ufgZzsnJvsnOvgnKsHPAq0xPfoWTBdgW/06xx8f5+09G/P/5hes0MvzKwGvlGvTf75QPfi+7Os45yrDew+4lwRESlFIVkIFWI0cJtzLh74G/8b+XkMGG5mKcBU4DZv4omIeMc5txt4BN98mEvMrJqZVTKz88zs//kPexw43cyeNLO6ZhZlZrfhewzuviMu+Q5wK9CDP8+TKZB/RKWVmb0M9PLfr6hjI4GKQAV/6+5K+fbH+49pgO///5/7R4qCKQpfMbkD32N4BRVt55vZGf5H9v4BzHHObfCfmw2kAhXN7BH+9490IiLigbAohPz/6nY6vomni4D/8L/HOy4HxjnnmgLnA++aWVh8bhGRo+Gfl3IXvkeGU/E9xnUrMMW/fxVwBr65Kcn4Hju7FOjnnJt1xOU+xjfa8Z1zrqB5Rfl1N7N0YA/wI74CoItzbkkR54zA12b7dXxNFTKAMfn2vwjsAn7D98je9cVkOB7O/9/xwDpgI75H8eYUcOz7+Eao0oB4fKNmAN8AXwO/+69xgHyP0YmISOkz51zxR5VBZtYS+NI518E/Wfc359yfnm03s2VAf/+/yGFmfwDdnHPbSjOviIiEHjP7BJjhnHshgGPHASnOuYeCHkxERI5bWIyMOOf2AGvNbAiA+RzqtrMe6OPf3g5f+9hUT4KKiEjIMLMm+EbIEr3OIiIiJS8kCyEz+wDfxN+2/oXxrsW3dsO1ZvYrsAz/SunA3cD1/u0fACNdqA6DiYhIqTCzW4CFwNvOuZ+9ziMiIiUvZB+NExEREREROVYhOSIkIiIiIiJyPEJuobb69eu7li1beh1DRERERETKqKSkpO3OuQZFHRNyhVDLli1JTNS8VRERERERKZiZrSvuGD0aJyIiIiIi5Y4KIRERERERKXdUCImIiIiISLkTcnOECpKVlUVKSgoHDhzwOoqUQZGRkTRt2pRKlSp5HUVEREREyoiwKIRSUlKIioqiZcuWmJnXcaQMcc6xY8cOUlJSaNWqlddxRERERKSMCItH4w4cOEC9evVUBMmfmBn16tXTaKGIiIiIHCYsCiFARZAUSj8bIiIiInKksCmEREREREREAqVCqIRs2bKFyy67jNatWxMfH8/555/P77//HtR79urVq9jFZV944QX279+f9/78889n165dx33vli1bsn379iKPWblyJXFxcXTq1Ik1a9Yc1fV//PFHZs+eHZRcIiIiIiIqhEqAc46BAwfSq1cv1qxZQ1JSEk899RRbt271OtqfCqGpU6dSu3btUrn3lClTGDx4MAsXLqR169ZHde6xFkIiIiIiIoFQIVQCfvjhBypVqsRNN92Uty02NpaePXvy448/cuGFF+Ztv/XWWxk3bhzgG70YNWoUcXFxJCQksGDBAvr160fr1q154403AIo8P7+bb76ZhIQEoqOjefTRRwF46aWX2LRpE71796Z3795599y+fTv3338/r776at75jz32GM8++ywAzzzzDF26dCEmJibvWoVJTk6mXbt2XH/99URHR9O3b18yMjKYOnUqL7zwAq+//nrevSdMmEDXrl2Ji4vjxhtvJCcnB4Cvv/6azp07ExsbS58+fUhOTuaNN97g+eefJy4ujpkzZ5Kamsqll15Kly5d6NKlC7NmzQJgx44d9O3bl+joaK677jqcc8V/w0RERESk3AuL9tn53TFlKYs27SnRa8Y1rskLl3QodP/SpUuJj48/pms3b96cRYsWceeddzJy5EhmzZrFgQMH6NChw2GFVXGefPJJ6tatS05ODn369GHx4sX83//9H8899xw//PAD9evXP+z4YcOGcccdd/DXv/4VgEmTJvHNN98wffp0Vq1axbx583DOMWDAAGbMmMGZZ55Z6L1XrVrFBx98wJgxYxg6dCgff/wxw4cP56abbqJGjRr87W9/Y8WKFUycOJFZs2ZRqVIlbrnlFt577z3OO+88rr/+embMmEGrVq1IS0ujbt26h50LcMUVV3DnnXdyxhlnsH79evr168eKFSt4/PHHOeOMM3jkkUf46quveOutt47huyAiIiIi5U3YFUKhZsCAAQB07NiR9PR0oqKiiIqKokqVKkc1l2fSpEmMHj2a7OxsNm/ezPLly4mJiSn0+E6dOrFt2zY2bdpEamoqderUoVmzZrz44otMnz6dTp06AZCens6qVauKLIRatWpFXFwcAPHx8SQnJ//pmO+++46kpCS6dOkCQEZGBieccAJz5szhzDPPzFvjp27dugXe49tvv2X58uV57/fs2UN6ejozZszgk08+AeCCCy6gTp06RfwpiYiIiIj4hF0hVNTITbBER0czefLkAvdVrFiR3NzcvPdHrmdTpUoVACpUqJD3+tD77OzsYs8HWLt2Lc8++yzz58+nTp06jBw5MqB1c4YMGcLkyZPZsmULw4YNA3zznUaNGsWNN95Y7PlHfgaAiIgIMjIy/nSMc46rr76ap5566rDtX3zxRUD3yM3NZc6cOURGRgacS0RERESkMEGbI2RmY81sm5ktLea4LmaWbWaDg5Ul2M4++2wyMzMZPXp03rbFixczc+ZMWrRowfLly8nMzGTXrl189913R3XtQM7fs2cP1atXp1atWmzdupVp06bl7YuKimLv3r0FXnvYsGF8+OGHTJ48mSFDhgDQr18/xo4dS3p6OgAbN25k27ZtR5W5IH369GHy5Ml510pLS2PdunV069aNGTNmsHbt2rztBeXu27cvL7/8ct77RYsWAXDmmWfy/vvvAzBt2jR27tx53FlFRERE5Ogs2rib9Mxsr2MclWA2SxgH9C/qADOLAP4FTA9ijqAzMz799FO+/fZbWrduTXR0NKNGjaJhw4Y0a9aMoUOH0qFDB4YOHZr3yFmgAjk/NjaWTp06ceqpp3LFFVfQo0ePvH033HAD/fv3z2tYkF90dDR79+6lSZMmNGrUCPAVHFdccQXdu3enY8eODB48uNBC6mi0b9+eJ554gr59+xITE8O5557L5s2badCgAaNHj2bQoEHExsbmjUxddNFFfPrpp3nNEl566SUSExOJiYmhffv2ec0kHn30UWbMmEF0dDSffPIJzZs3P+6sIiIiIhK4nFzHxW/P54oJC7yOclQsmF22zKwl8KVzrsDn1czsDiAL6OI/ruDny/JJSEhwR66ds2LFCtq1a3fceSV86WdEREREJDi+WbmN/mPmMnFEPEPjGnsdBwAzS3LOJRR1jGfts82sCTAQeD2AY28ws0QzS0xNTQ1+OBERERERCchb89ZTr1olLu5wotdRjoqX6wi9ANznnMst7kDn3GjnXIJzLqFBgwalEE1ERERERIqTmp7JlKVbuCqhGVUqRngd56h42TUuAfjQzADqA+ebWbZzboqHmUREREREJEDvJqWQleO49rTQm6ftWSHknGt16LWZjcM3R0hFkIiIiIhICHDO8dbc9XRrUYfohlFexzlqQSuEzOwDoBdQ38xSgEeBSgDOuTeCdV8REREREQm+Oet2snxrOm8OjfU6yjEJWiHknLv8KI4dGawcIiIiIiJS8t6cu54aVSIYVkY6xR0tL5slhJWUlBQuvvhi2rRpQ+vWrbn99ts5ePAgAOPGjePWW2/1OOGf1ahR46i25zdz5kyio6OJi4sjIyPjqO47ZcoUli9fflTnBJpLRERERIJv74FsJi7axLDYJtSo4mXbgWOnQqgEOOcYNGgQl1xyCatWreL3338nPT2dBx98MGj3zM72duXe9957j1GjRrFo0SKqVq16VOceayEkIiIiImXDxEUb2Xcwh+u6hV6ThENUCJWA77//nsjISK655hoAIiIieP755xk7diz79+8HYMOGDfTq1Ys2bdrw+OOPA7Bv3z4uuOACYmNj6dChAxMnTgQgKSmJs846i/j4ePr168fmzZsB6NWrF3fccQcJCQk8+eSTtGjRgtzc3LxrNWvWjKysLNasWUP//v2Jj4+nZ8+erFy5EoC1a9fSvXt3OnbsyEMPPVTs5/rxxx/p1asXgwcP5tRTT+XKK6/EOcebb77JpEmTePjhh7nyyisBeOaZZ+jSpQsxMTE8+uijedcYP348MTExxMbGMmLECGbPns3nn3/OPffcQ1xcHGvWrCmxvCIiIiJSOt6cu57ohlGc1ry211GOWWiOYxXhjrmfsShtY4leM65uE1447eJC9y9btoz4+PjDttWsWZPmzZuzevVqAObNm8fSpUupVq0aXbp04YILLmDdunU0btyYr776CoDdu3eTlZXFbbfdxmeffUaDBg2YOHEiDz74IGPHjgXg4MGDJCYmArBgwQJ++uknevfuzZdffkm/fv2oVKkSN9xwA2+88QZt2rRh7ty53HLLLXz//ffcfvvt3HzzzVx11VW8+uqrAX32hQsXsmzZMho3bkyPHj2YNWsW1113HT///DMXXnghgwcPZvr06axatYp58+bhnGPAgAHMmDGDevXq8cQTTzB79mzq169PWloadevWZcCAAXnnAvTp06fE8oqIiIhIcC3dvIe563fx3ID2+JfCCUlhVwiVVeeeey716tUDYNCgQfz888+cf/753H333dx3331ceOGF9OzZk6VLl7J06VLOPfdcAHJycmjUqFHedYYNG3bY64kTJ9K7d28+/PBDbrnlFtLT05k9ezZDhgzJOy4zMxOAWbNm8fHHHwMwYsQI7rvvvmJzd+3alaZNmwIQFxdHcnIyZ5xxxmHHTJ8+nenTp9OpUycA0tPTWbVqFb/++itDhgyhfv36ANStW/dP1y/pvCIiIiISXG/NW0+lCGNEfFOvoxyXsCuEihq5CZb27dszefLkw7bt2bOH9evXc/LJJ7NgwYI/VctmximnnMKCBQuYOnUqDz30EH369GHgwIFER0fzyy+/FHiv6tWr570eMGAADzzwAGlpaSQlJXH22Wezb98+ateuzaJFiwo8/2ir9ipVquS9joiIKHBuknOOUaNGceONNx62/eWXXy72+rm5uSWaV0RERESCJzM7h/GJKQzs0Ij6NaoUf0IZpjlCJaBPnz7s37+f8ePHA75RnLvvvpuRI0dSrVo1AP773/+SlpZGRkYGU6ZMoUePHmzatIlq1aoxfPhw7rnnHhYsWEDbtm1JTU3NK4SysrJYtmxZgfetUaMGXbp04fbbb+fCCy8kIiKCmjVr0qpVKz766CPAV6T8+uuvAPTo0YMPP/wQ8DU7KCn9+vVj7NixpKenA7Bx40a2bdvG2WefzUcffcSOHTsASEtLAyAqKoq9e/cCeJJXRERERI7NlCVbSNufxXWnhW6ThENUCJUAM+PTTz/lo48+ok2bNpxyyilERkbyz3/+M++Yrl27cumllxITE8Oll15KQkICS5YsoWvXrsTFxfH444/z0EMPUblyZSZPnsx9991HbGwscXFxzJ49u9B7Dxs2jAkTJhz2yNx7773HW2+9RWxsLNHR0Xz22WcAvPjii7z66qt07NiRjRtLbh5V3759ueKKK/IaGwwePJi9e/cSHR3Ngw8+yFlnnUVsbCx33XUXAJdddhnPPPMMnTp1Ys2aNaWeV0RERESOzZtz19OiTlX6tKnvdZTjZs45rzMclYSEBHeoWcAhK1asoF27dh4lklCgnxERERGR47N2x35O+ud3PN6vLY/0PcXrOEUysyTnXEJRx2hESEREREREivX2/PWYwTVdmnkdpUSoEBIRERERkSLl5DrenreBfm0b0KxOVa/jlAgVQiIiIiIiUqTpv20jZfeBsGiScIgKIRERERERKdKbc9fToEZlLmrf0OsoJUaFkIiIiIiIFGrr3kw+X7aVqxOaUbli+JQP4fNJRERERESkxI1P3EB2ruParuHRJOEQFUIlJCIigri4ODp06MCQIUPYv39/ocfu2rWL11577bjuN27cODZt2pT3/rrrrmP58uVFntOrVy+ObD1e1Pb8UlNTOe200+jUqRMzZ848qqyLFi1i6tSpR3VOoLlEREREJHicc4ydt4HTW9bh1BOjvI5TolQIlZCqVauyaNEili5dSuXKlXnjjTcKPTYYhdCbb75J+/btj+uaRfnuu+/o2LEjCxcupGfPnkd17rEWQiIiIiLirfkbdrFyWzojw6Rldn4qhIKgZ8+erF69GoDnnnuODh060KFDB1544QUA7r//ftasWUNcXBz33HMPAM888wxdunQhJiaGRx99FIDk5GTatWvH9ddfT3R0NH379iUjI4PJkyeTmJjIlVdeSVxcHBkZGYeNntx8880kJCQQHR2dd61A1ahRgwcffJDY2Fi6devG1q1bWbRoEffeey+fffZZ3v2mT59O9+7d6dy5M0OGDCE9PR2A+fPnc/rppxMbG0vXrl3ZvXs3jzzyCBMnTiQuLo6JEyeyb98+/vKXv9C1a1c6derEZ599BkBGRgaXXXYZ7dq1Y+DAgWRkZBz/N0NEREREjtn4xBQiK1ZgaGxjr6OUuIpeByhpW967gwPrF5XoNSObx9HwyhcCOjY7O5tp06bRv39/kpKSePvtt5k7dy7OOU477TTOOussnn76aZYuXcqiRb6c06dPZ9WqVcybNw/nHAMGDGDGjBk0b96cVatW8cEHHzBmzBiGDh3Kxx9/zPDhw3nllVd49tlnSUj484K5Tz75JHXr1iUnJ4c+ffqwePFiYmJiAsq/b98+unXrxpNPPsm9997LmDFjeOihh/j73/9OYmIir7zyCtu3b+eJJ57g22+/pXr16vzrX//iueee4/7772fYsGFMnDiRLl26sGfPHqpVq3bYuQAPPPAAZ599NmPHjmXXrl107dqVc845h//85z9Uq1aNFStWsHjxYjp37hzgd0hERERESlpmdg4fLNzIJR0aUqtqJa/jlLiwK4S8kpGRQVxcHOAbEbr22mt5/fXXGThwINWrVwdg0KBBzJw5kwEDBhx27vTp05k+fTqdOnUCID09nVWrVtG8eXNatWqVd934+HiSk5OLzTJp0iRGjx5NdnY2mzdvZvny5QEXQpUrV+bCCy/Mu99///vfPx0zZ84cli9fTo8ePQA4ePAg3bt357fffqNRo0Z06dIFgJo1axZ4j+nTp/P555/z7LPPAnDgwAHWr1/PjBkz+L//+z8AYmJiAs4sIiIiIiXvq+XbSNufxdVh+FgchGEhFOjITUk7NEfoWDjnGDVqFDfeeONh25OTk6lSpUre+4iIiGIfF1u7di3PPvss8+fPp06dOowcOZIDBw4EnKVSpUqYWd79srOzC8x77rnn8sEHHxy2fcmSJQHdwznHxx9/TNu2bQPOJSIiIiKl653EDTSqWYVzT2ngdZSg0ByhIOrZsydTpkxh//797Nu3j08//ZSePXsSFRXF3r17847r168fY8eOzZtns3HjRrZt21bktY+8xiF79uyhevXq1KpVi61btzJt2rSS/VBAt27dmDVrVt48qH379vH777/Ttm1bNm/ezPz58wHYu3cv2dnZBX7el19+GeccAAsXLgTgzDPP5P333wdg6dKlLF68uMSzi4iIiEjxUtMzmbpiG8M7NyWignkdJyjCbkSoLOncuTMjR46ka9eugK/F9aHH33r06EGHDh0477zzeOaZZ1ixYgXdu3cHfA0LJkyYQERERKHXHjlyJDfddBNVq1bll19+ydseGxtLp06dOPXUU2nWrFne42slqUGDBowbN47LL7+czMxMAJ544glOOeUUJk6cyG233UZGRgZVq1bl22+/pXfv3jz99NPExcUxatQoHn74Ye644w5iYmLIzc2lVatWfPnll9x8881cc801tGvXjnbt2hEfH1/i2UVERESkeO8v2Eh2ruOqhKZeRwkaO/Sv8qEiISHBHbm2zIoVK2jXrp1HiSQU6GdEREREJHDxz88AIOnOMz1OcmzMLMk59+euYvno0TgREREREcmzdPMeFqTs5uowHg0CFUIiIiIiIpLPO4kpVKxgXN6piddRgipsCqFQe8RPSo9+NkREREQCk52Ty4SkFC5odwINalQp/oQQFhaFUGRkJDt27NAvvPInzjl27NhBZGSk11FEREREyrz//p7Klr2ZYbt2UH5h0TWuadOmpKSkkJqa6nUUKYMiIyNp2jS8n3EVERERKQnvJKZQt1olLmh3otdRgi4sCqFKlSrRqlUrr2OIiIiIiISsXRlZTFm6hetOa07limHx4FiRwv8TioiIiIhIsT76dROZ2blcnRD+j8WBCiEREREREQHemb+BdifWIKFZLa+jlAoVQiIiIiIi5dzq7fuYlbyTqxOaYWZexykVKoRERERERMq58YkbqGAwPD681w7KT4WQiIiIiEg5lpvrGJ+YwjltGtCkVlWv45QaFUIiIiIiIuXYzLU7WLczg6u7lK/lRlQIiYiIiIiUY+/MTyGqSkUu6dDQ6yilSoWQiIiIiEg5tS8zm48Wb2JIbCOqVQ6LJUYDpkJIRERERKSc+nTpFtIzc8rN2kH5qRASERERESmnxiduoFXdapzRqq7XUUpd0AohMxtrZtvMbGkh+680s8VmtsTMZptZbLCyiIiIiIjI4TbtPsB3q7YzPL4JFSqUj7WD8gvmiNA4oH8R+9cCZznnOgL/AEYHMYuIiIiIiOTzwcKN5DoYHl++usUdErQZUc65GWbWsoj9s/O9nQOUz++AiIiIiIgHJiSl0LV5bU5pUMPrKJ4oK3OErgWmFbbTzG4ws0QzS0xNTS3FWCIiIiIi4Wfp5j0s2rSH4Z3L71iE54WQmfXGVwjdV9gxzrnRzrkE51xCgwYNSi+ciIiIiEgYmpC0kYgKxmWdGnsdxTOeNgs3sxjgTeA859wOL7OIiIiIiJQHubmO9xak0L9tAxrUqOJ1HM94NiJkZs2BT4ARzrnfvcohIiIiIlKe/LhmBym7DzCinDZJOCRoI0Jm9gHQC6hvZinAo0AlAOfcG8AjQD3gNTMDyHbOJQQrj4iIiIiI+JokRFWpyIAODb2O4qlgdo27vJj91wHXBev+IiIiIiJyuP0Hs5m8eDODYxpRtVKE13E85XmzBBERERERKR1fLNvK3sxsRiSU78fiQIWQiIiIiEi58W5SCk1rRXLWSfW8juI5FUIiIiIiIuVAanomX/+WypWdm1Khgnkdx3MqhEREREREyoEPF24iJ9cxPL6J11HKBBVCIiIiIiLlwLtJKcQ1rkmHRjW9jlImqBASEREREQlzv21LZ/6GXQwv52sH5adCSEREREQkzE1ISqGCweWd9FjcISqERERERETCmHOOCQtS6NOmPo1rRXodp8xQISQiIiIiEsZmrU0jOS2DEXos7jAqhEREREREwtiEBRupVjmCgR0beR2lTFEhJCIiIiISpjKzc5i0aBMDOzSkRpWKXscpU1QIiYiIiIiEqa+Wb2NnRpa6xRVAhZCIiIiISJiasCCFE6OqcE6b+l5HKXNUCImIiIiIhKG0/Qf5cvHz0IQAACAASURBVPlWLu/UmIoR+rX/SPoTEREREREJQx/9uomsHKducYVQISQiIiIiEobGJ6bQ7sQadGpSy+soZZIKIRERERGRMLN6+z5mJ+/k6oRmmJnXccokFUIiIiIiImFmQlIKZnBl5yZeRymzVAiJiIiIiIQR5xzjE1Poc3J9mtau6nWcMkuFkIiIiIhIGJm1No21afu5KkFNEoqiQkhEREREJIyMT0qheuUIBnZs5HWUMk2FkIiIiIhImMjIymHSok1cGtOIGlUqeh2nTFMhJCIiIiISJr5YtpXdB7K1dlAAVAiJiIiIiISJ8YkbaFIrkt4n1/c6SpmnQkhEREREJAxs3ZvJ17+lMrxzUyIqaO2g4qgQEhEREREJAx8s3EhOrmOEusUFRIWQiIiIiEgYeDcphfimtYhuGOV1lJCgQkhEREREJMQt3byHBSm7tXbQUVAhJCIiIiIS4t5NSqFiBeOyuCZeRwkZKoREREREREJYTq5jQtJGzjv1BE6IquJ1nJChQkhEREREJIR9v2o7m/Yc0GNxR0mFkIiIiIhICBuftIHaVStxYfsTvY4SUlQIiYiIiIiEqL0HsvlkyRaGxjYislKE13FCigohEREREZEQ9cmSzew/mMNVCc28jhJyVAiJiIiIiISo8YkpnFSvGqe3rON1lJCjQkhEREREJARt2JnBD2u2c1V8U8zM6zghR4WQiIiIiEgIem9BCs7BCHWLOyYqhEREREREQoxzjvFJKZzRqi4n1avudZyQpEJIRERERCTEJKXsZsXWdK0ddBxUCImIiIiIhJh35m+gSsUKDIlt7HWUkBW0QsjMxprZNjNbWsh+M7OXzGy1mS02s87ByiIiIiIiEi4ys3N4f+FGBnZoSO2qlbyOE7KCOSI0DuhfxP7zgDb+rxuA14OYRUREREQkLHy5fCtp+7MY2UVrBx2PoBVCzrkZQFoRh1wMjHc+c4DaZtYoWHlERERERMLBuPkpNK4ZyTmnNPA6Skjzco5QE2BDvvcp/m1/YmY3mFmimSWmpqaWSjgRERERkbJmy54DTFu5jasSmhJRQWsHHY+QaJbgnBvtnEtwziU0aKDKV0RERETKp/cWbCQn13G1usUdNy8LoY1A/gcbm/q3iYiIiIjIEZxzjJu/gW4t6nDqiVFexwl5XhZCnwNX+bvHdQN2O+c2e5hHRERERKTMWpCym6Vb9jKyi0aDSkLFYF3YzD4AegH1zSwFeBSoBOCcewOYCpwPrAb2A9cEK4uIiIiISKgb5187aFhcgdPq5SgFrRByzl1ezH4H/DVY9xcRERERCRdaO6jkhUSzBBERERGR8kxrB5U8FUIiIiIiImWc1g4qeSqERERERETKMK0dFBwqhEREREREyjCtHRQcKoRERERERMoorR0UPCqERERERETKKK0dFDwqhEREREREyiitHRQ8KoRERERERMogrR0UXCqERERERETKoK+Wb9PaQUFUbCFkZrebWU3zecvMFphZ39IIJyIiIiJSXo2bv0FrBwVRICNCf3HO7QH6AnWAEcDTQU0lIiIiIlKObd2bydSV2xgRr7WDgiWQQujQn/z5wLvOuWX5tomIiIiISAl7b0GKb+0gdYsLmkAKoSQzm46vEPrGzKKA3ODGEhEREREpn5xzvD1vA6c1r007rR0UNBUDOOZaIA74wzm338zqAdcEN5aIiIiISPl0aO2g1y/t6HWUsBbIiNB/nXMLnHO7AJxzO4DngxtLRERERKR8GjN3PVUrVeCyTlo7KJgKHREys0igGlDfzOrwv3lBNQF9V0RERERESlh6ZjbvLUhhaGxjrR0UZEU9GncjcAfQGEjif4XQHuCVIOcSERERESl3Ji7aRHpmDjd0a+F1lLBXaCHknHsReNHMbnPOvVyKmUREREREyqXRc9YR3TCK7i3reB0l7BXbLME597KZnQ60zH+8c258EHOJiIiIiJQrv27azbz1u3jh4mjMtFpNsBVbCJnZu0BrYBGQ49/sABVCIiIiIiIlZMyc9VSpWIERCVo7qDQE0j47AWjvnHPBDiMiIiIiUh7tP5jNu0kpDI5pRN1qlb2OUy4E0j57KdAw2EFERERERMqrSYs2s+dAtpoklKJARoTqA8vNbB6QeWijc25A0FKJiIiIiJQjY+auo22D6vQ8qa7XUcqNQAqhx4IdQkRERESkvFq6eQ+zk3fy7EXt1SShFAXSNe4nM2sBtHHOfWtm1YCI4EcTEREREQl/Y+aup3JEBa5Wk4RSVewcITO7HpgM/Me/qQkwJZihRERERETKg4ysHN5NTGFQx4bUr1HF6zjlSiDNEv4K9AD2ADjnVgEnBDOUiIiIiEh58PHizezMyOJ6NUkodYEUQpnOuYOH3phZRXzrCImIiIiIyHEYM2cdJ9evTq/W9byOUu4EUgj9ZGYPAFXN7FzgI+CL4MYSEREREQlvK7fuZcYfaVx3WnMqVFCThNIWSCF0P5AKLAFuBKYCDwUzlIiIiIhIuBszdz0VKxgjuzTzOkq5FEjXuFxgjP9LRERERESOU2Z2Du/M38AlHRpyYpSaJHih0ELIzCY554aa2RIKmBPknIsJajIRERERkTD16ZIt7NifxfXdmnsdpdwqakTodv9/LyyNICIiIiIi5cXoOetoWbcq57Rp4HWUcqvQOULOuc3+l7c459bl/wJuKZ14IiIiIiLhZVVqOj+s3sH1p7VQkwQPBdIs4dwCtp1X0kFERERERMqDN+euJ6KCcU1XNUnwUlFzhG7GN/JzkpktzrcrCpgV7GAiIiIiIuHmYHYub8/fwEXtT6RRzUiv45RrRc0Reh+YBjyFr4X2IXudc2lBTSUiIiIiEoY+XbKZ1PSD3KAmCZ4rtBByzu0GdgOXA5jZCUAkUMPMajjn1pdORBERERGR8PDa7GRa1a1Gv7YneB2l3Ct2jpCZXWRmq4C1wE9AMr6RIhERERERCdDSzXuY8UcaN5+uJgllQSDNEp4AugG/O+daAX2AOUFNJSIiIiISZl6fvY4qFStwTRc1SSgLAimEspxzO4AKZlbBOfcDkBDIxc2sv5n9Zmarzez+AvY3N7MfzGyhmS02s/OPMr+IiIiISJm390A245M2MCyuMfVrVPE6jlB0s4RDdplZDWAm8J6ZbQP2FXeSmUUAr+Jrv50CzDezz51zy/Md9hAwyTn3upm1B6YCLY/yM4iIiIiIlGkTFqSQnpnDLae39DqK+AUyInQxkAHcAXwNrAEuCuC8rsBq59wfzrmDwIf+a+XngJr+17WATYGEFhEREREJFc45XpuVTOemtejavLbXccSv2ELIObcPaACcD6ThG8HZEcC1mwAb8r1P8W/L7zFguJml4BsNuq2gC5nZDWaWaGaJqampAdxaRERERKRs+HltGku37OWW01tipiYJZUUgXeOuA+YBg4DBwBwz+0sJ3f9yYJxzrim+QutdM/tTJufcaOdcgnMuoUGDBiV0axERERGR4HttVjK1IityeafGXkeRfAKZI3QP0OnQKJCZ1QNmA2OLOW8jkL8lRlP/tvyuBfoDOOd+MbNIoD6wLYBcIiIiIiJl2pY9B/h4yWb+2qMl1SoH8qu3lJZA5gjtAPbme7/Xv60484E2ZtbKzCoDlwGfH3HMenztuDGzdvgWbNWzbyIiIiISFt6at56sHMdN3Vt6HUWOEEhZuhqYa2af4WtucDGw2MzuAnDOPVfQSc65bDO7FfgGiADGOueWmdnfgUTn3OfA3cAYM7vTf+2Rzjl33J9KRERERMRj2Tm5/OeXdZzTpj5tT6jhdRw5QiCF0Br/1yGf+f8bVdyJzrmp+Jog5N/2SL7Xy4EeAWQQEREREQkpX63YxoZdB3jxkg5eR5ECFFsIOeceL40gIiIiIiLh5LVZyTSpFclF7U/0OooUIJA5QiIiIiIichRWpaYz/fdUbuzegooR+pW7LNJ3RURERESkhL3xyzoqVjCuO62511GkECqERERERERK0P6D2bw9bwODOjaiUc1Ir+NIIQqdI2RmFfGt8zMQOLT600Z8zRLecs5lBT+eiIiIiEhombhoEzszsrilRwuvo0gRimqW8C6wC3gMSPFvawpcDUwAhgU1mYiIiIhICHptdjLRDaM486R6XkeRIhRVCMU75045YlsKMMfMfg9iJhERERGRkDR//S4SN+zm1UEdMTOv40gRipojlGZmQ8ws7xgzq2Bmw4CdwY8mIiIiIhJaXpudTI0qEQyPb+J1FClGUYXQZcBgYKuZ/W5mq4AtwCD/PhERERER8UtNz+TDhRsZ3rkpNSMreR1HilHoo3HOuWT884DMrJ5/247SiSUiIiIiElpem5XMgexcbu/ZyusoEoCiusYNKmBb3mvn3CdByiQiIiIiElIysnJ4ZVYyF7U/kVNPjPI6jgSgqGYJFxWxzwEqhEREREREgPGJG9i+7yB/69Xa6ygSoKIKoTHAL845V1phRERERERCTU6u498//kGXZrXpeVJdr+NIgIpqljACSDKzD81spJk1LK1QIiIiIiKh4otlW1i1fR/39G6tltkhpKhmCTcDmNmpwHnAODOrBfwAfA3Mcs7llEpKEREREZEy6tkf19CqbjUGdtC4QSgpakQIAOfcSufc8865/sDZwM/AEGBusMOJiIiIiJRlvySnMSt5J3eeeRIVI4r91VrKkKK6xhX2gOMc4Cfn3L7gRBIRERERCQ3//ukP6lStxDVdm3kdRY5SUc0SkgrZ7oCK5nsA8n7n3HslH0tEREREpGxbvX0fnyzZzKizT6ZGlaJ+rZayqKjvWFvn3MHCdppZA+AnQIWQiIiIiJQ7z//0B5UqVOC2M7SAaigqqhCabWYp+BojfO2cS86/0zmXamb3BTOciIiIiEhZtD09k7fnr2dEfFMa1oz0Oo4cg6K6xiWYWUugP/CCmTXB1yhhGr45QpnOuS9KJaWIiIiISBny+i/ryMjK5a6zTvI6ihyjIltbOOeSnXNvOOcuAU4HvgDOAWaa2VelEVBEREREpCw5kJXDyz+v5YJ2J9C+YZTXceQYBTyryzmXBXzv/8I/QiQiIiIiUq68m5RCavpB/tartddR5DgU1T57Cb4OcQVyzsUEJZGIiIiISBmVm+v4949riG9ai7Na1/M6jhyHokaELiy1FCIiIiIiIeDL5Vv5LXUfHw7vjG81GQlVRTVLWFeaQUREREREyrpnf1pDizpVuTSmkddR5DgV2SxBRERERER85q7bycw/0rjzzJOoGKFfo0OdvoMiIiIiIgF49sc11K5aib90be51FCkBhRZCZvad/7//Kr04IiIiIiJlz8qte/l4yWZuPr0FUZEBN16WMqyo72IjMzsdGGBmHwKHzQZzzi0IajIRERERkTLin9+tpmqlCO48UwuohouiCqFHgIeBpsBzR+xzwNnBCiUiIiIiUlas2b6P9xdu5I6erWhQo4rXcaSEFNU1bjIw2cweds79oxQziYiIiIiUGU99t5pKFUwLqIaZYh9wdM79w8wGAGf6N/3onPsyuLFERERERLyXnLafdxI3cPPpLWlYM9LrOFKCiu0aZ2ZPAbcDy/1ft5vZP4MdTERERETEa//6fjUVzLi3t0aDwk0gLS8uAOKcc7kAZvYOsBB4IJjBRERERES8tHF3BmPnbeCars1oWruq13GkhAW6jlDtfK9rBSOIiIiIiEhZ8v9+WEOuc9x/9sleR5EgCGRE6ClgoZn9gK+F9pnA/UFNJSIiIiLioS17DjD6l3WMiG9Ky7rVvI4jQRBIs4QPzOxHoIt/033OuS1BTSUiIiIi4qF///QHB3NyeeCcNl5HkSAJaFlc59xm4PMgZxERERER8VxqeiavzU7m8k5NOLl+da/jSJAEOkdIRERERKRceH7GH2Rk5fCgRoPCmgohERERERG/tP0HeeXnZAbHNKLdiVFex5EgKrIQMrMIM1t5rBc3s/5m9puZrTazAhssmNlQM1tuZsvM7P1jvZeIiIiIyPF6aeZa9mZm89A5p3gdRYKsyDlCzrkcfyHT3Dm3/mgubGYRwKvAuUAKMN/MPnfOLc93TBtgFNDDObfTzE44+o8gIiIiInL89hzI4sWZa7mkQ0NiGtf0Oo4EWSDNEuoAy8xsHrDv0Ebn3IBizusKrHbO/QFgZh8CFwPL8x1zPfCqc26n/5rbjiK7iIiIiEiJeeXnZHZlZPGQ5gaVC4EUQg8f47WbABvyvU8BTjvimFMAzGwWEAE85pz7+sgLmdkNwA0AzZs3P8Y4IiIiIiIFS8/M5rmf1nB+uxOIb1bb6zhSCgJZR+gnM2sBtHHOfWtm1fAVLSV1/zZAL6ApMMPMOjrndh2RYTQwGiAhIcGV0L1FRERERAB4fXYyO/Zn8fC5mhtUXhTbNc7MrgcmA//xb2oCTAng2huBZvneN/Vvyy8F+Nw5l+WcWwv8jq8wEhEREREpFbszsvjX96s595T6dGtRx+s4UkoCaZ/9V6AHsAfAObcKCKSpwXygjZm1MrPKwGX8eVHWKfhGgzCz+vgelfsjoOQiIiIiIiXgmR/XsGN/Fv88v53XUaQUBVIIZTrnDh56Y2YVgWIfT3POZQO3At8AK4BJzrllZvZ3MzvUaOEbYIeZLQd+AO5xzu042g8hIiIiInIsNu85wHM/rWFYXGMSNDeoXAmkWcJPZvYAUNXMzgVuAb4I5OLOuanA1CO2PZLvtQPu8n+JiIiIiJSqx775jexcx5Pnnep1FCllgYwI3Q+kAkuAG/EVNg8FM5SIiIiISLCt3LqXt+Zt4KbuLWldv7rXcaSUBdI1LtfM3gHm4nsk7jf/SI6IiIiISMh6YNpKqlaqoHWDyqlAusZdAKwBXgJeAVab2XnBDiYiIiIiEiyz16bx6ZIt3Nv7ZE6IquJ1HPFAIHOE/g30ds6tBjCz1sBXwLRgBhMRERERCQbnHPd+uZwTo6pw15kneR1HPBJIIbT3UBHk9wewN0h5RERERESC6otlW5mVvJPXL+1I9SqB/Dos4ajQ77yZDfK/TDSzqcAkfHOEhuBbI0hEREREJKRk5+Ry/9QVnNKgOtee1tzrOOKhokrgi/K93gqc5X+dClQNWiIRERERkSB5JzGFFVvT+fjqBCpFBNJAWcJVoYWQc+6a0gwiIiIiIhJM+w9m88jXv9GtRR0GdmzodRzxWLEPRZpZK+A2oGX+451zA4IXS0RERESkZL04cy2b9hzgg+GdMTOv44jHApkdNgV4C/gCyA1uHBERERGRkrdj30Ge/n41F7Y/kTNb1/M6jpQBgRRCB5xzLwU9iYiIiIhIkDz57SrSM7N5+oJ2XkeRMiKQQuhFM3sUmA5kHtronFsQtFQiIiIiIiUkOW0/r85KZmSXZkQ3jPI6jpQRgRRCHYERwNn879E4538vIiIiIlKmPTB1JRUMHu/X1usoUoYEUggNAU5yzh0MdhgRERERkZL08x87+GDhRh46pw1Na2sFGPmfQJqnLwVqBzuIiIiIiEhJysl13PbpUprVjmRUn5O9jiNlTCAjQrWBlWY2n8PnCKl9toiIiIiUWWPmrGPRpj1MHBFPtcqB/Nor5UkgPxGPBj2FiIiIiEgJStt/kAenreSs1vUYEtvI6zhSBhVbCDnnfiqNICIiIiIiJeWRr39jV0YWL13SQYunSoGKLYTMbC++LnEAlYFKwD7nXM1gBhMRERERORaLN+3h9dnJ3Hx6S2Ia61dWKVggI0J5zdbNV05fDHQLZigRERERkWPhnOO2T5dQp2ol/t5f7bKlcIF0jcvjfKYA/YKUR0RERETkmE1atIkZf6Tx5PmnUrdaZa/jSBkWyKNxg/K9rQAkAAeClkhERERE5Bjsy8zmb18sJ65xTa47rYXXcaSMC6Rr3EX5XmcDyfgejxMRERERKTOe/n41KbsP8MHwzkRUUIMEKVogc4SuKY0gIiIiIiLH6o8d+3jmxzVc2bkJZ5xUz+s4EgIKLYTM7JEiznPOuX8EIY+IiIiIyFG767NlVKxg/OvCdl5HkRBR1IjQvgK2VQeuBeoBKoRERERExHPfrNzGZ8u28tT5p9KkVlWv40iIKLQQcs79+9BrM4sCbgeuAT4E/l3YeSIiIiIipeVgdi63T1nKyfWrc+dZJ3kdR0JIkXOEzKwucBdwJfAO0Nk5t7M0gomIiIiIFOfFmX/wW+o+vry2K1UqRngdR0JIUXOEngEGAaOBjs659FJLJSIiIiJSjFWp6Tzy9W8MiD6RC9qf6HUcCTFFLah6N9AYeAjYZGZ7/F97zWxP6cQTEREREfmz3FzH9R8tpkrFCrx+aYzXcSQEFTVHqKgiSURERETEM6PnrOOnNTt4c2gsjWtFeh1HQpCKHREREREJKRt2ZnDvlyvo06Y+f+nazOs4EqJUCImIiIhIyHDOcePkxeQ4x5ghsZiZ15EkRBXZNU5EREREpCx5b8FGpq3cxgsXR9OqXjWv40gI04iQiIiIiISErXszuX3KUrq3qMOtZ7TyOo6EOBVCIiIiIhIS/u/TpaRn5vDWsFgiKuiRODk+KoREREREpMybsmQzk37dxCN929DuxCiv40gYUCEkIiIiImXazv0HufnjJcQ2rsm9vU/2Oo6ECTVLEBEREZEy7W9fLCd130G+uq4rlSL07/hSMvSTJCIiIiJl1n9/S2XsvA3c06s1nZvW9jqOhJGgFkJm1t/MfjOz1WZ2fxHHXWpmzswSgplHREREREJHemY213/0K20bVOfRvqd4HUfCTNAKITOLAF4FzgPaA5ebWfsCjosCbgfmBiuLiIiIiISeOz9bxvpdGbw5NJbIShFex5EwE8wRoa7AaufcH865g8CHwMUFHPcP4F/AgSBmEREREZEQMmnRJt6cu577zz6ZM06q53UcCUPBLISaABvyvU/xb8tjZp2BZs65r4q6kJndYGaJZpaYmppa8klFREREpMxITtvPDR/9SrcWdXi8X1uv40iY8qxZgplVAJ4D7i7uWOfcaOdcgnMuoUGDBsEPJyIiIiKeyMrJ5fIJC3DA+1d2Vpc4CZpg/mRtBJrle9/Uv+2QKKAD8KOZJQPdgM/VMEFERESk/Hrsm9+Ys24nowfH0KpeNa/jSBgLZiE0H2hjZq3MrDJwGfD5oZ3Oud3OufrOuZbOuZbAHGCAcy4xiJlEREREpIz6ftV2nvp+Ndd2bc6wTk2KP0HkOAStEHLOZQO3At8AK4BJzrllZvZ3MxsQrPuKiIiISOhJTc9k+PsLaNugBi9eEu11HCkHKgbz4s65qcDUI7Y9UsixvYKZRURERETKJucc13y4iLT9Wfz/9u47Po76zv/467u7WvUu2ZKrbGPccMPGxqYTIKYkQBJaQkkIISRHCneX3F2SHy0JB+GSCxdySQgX2gUCoYQeOj4wBmyDbbAN7kW2bMnqbet8f3/MSF4LyZaN5ZW07+fjsY9p35n5fHZH0nw0M999/htzyU7v01NUEaCPCyERERERkf35rzc28eyaan5z/lFMH5af7HAkRagbDhERERFJmvcrG/nhM2v4/JSh/MNxFckOR1KICiERERERSYqWcIyL/3cZpTlB/nTRDIwxyQ5JUohujRMRERGRpPjOEx+ybncrr14zj+LsYLLDkRSjK0IiIiIictjdt2Qb9y7Zxk9OG8/JR5QkOxxJQSqEREREROSwWrqtgW8+upJTjyjh+tOPTHY4kqJUCImIiIjIYVPdHOb8e5ZQlpvOw5cdTcCv01FJDj0jJCIiIiKHRTTucOEDy9jdGuGt7xxPSU56skOSFKZCSEREREQOix88vZqFG2p54MszmTlC3xckyaVrkSIiIiLS5x5Yuo073tjE908cw6WzRiQ7HBEVQiIiIiLSt5Zta+Dqv67k5HHF3H7O5GSHIwKoEBIRERGRPlTTEuYL9y2lNCfII5fPUucI0m/oGSERERER6ROxuMOF9y+jujnMm9ceR6k6RzgkbCxK++ZlhLetxMajYB2sdcBxsE4crNM5z5eeQ/akU0gfcRTGmGSH3q+oEBIRERGRPvHDZ9bw+oZa7r9kBrNGFiQ7nAHLiYZp3/gubR8tpO3jhbStewsbaTugbQTyy8iechrZU04ne8pppBUO66NoBw4VQiIiIiJyyP3vskr+8/828t0TxnDZ7JHJDmdAsfEYbWvfoHXN67R9vJD2DW9jo2EA0kdMpeDEK8mecCIZY+fgC2aBz4cxPvBexuf3hj5ijbtoXfUyLateouWDF2h863/d7Qyf4hZFR51O9sST8KVnJzPlpDDW2mTHcEBmz55tly5dmuwwRERERKQHb2ys5Yw/vM3c0YW89M1jSdNzQftlrSW0dTmNix6g8e0HiTfuAuMjY/RMsiacSPbEk8g88ngCOcUHvw/HIbxtJS2rXqL1w5doW/sGNhrCn1tK+RW/I++YLx7CjJLLGLPMWjt7n21UCImIiIjIofL2lnpO/8NiRuRnsvDb8xmSq+eC9iVat53GxX+m8a0HCFd+CP40cmecQ/78S8me/Bn8WX33fUtOpJ22j9+g+tEfEdq8jLxjL6Hsst98qmKrv+hNIaRb40RERETkkFi2rYEFd73N0Jx0XrlmnoqgHjihFpqWPUHjovtpXf0KWEvmEfMou/y/yZt74WErRHzBTHKmnkH2pFPY/ext1Dx5M21rXqP8a3eRO/NzhyWGZNIVIRERERH51FbuaOKU371FXkaAhd+ez6jCrGSH1O/EW+up/fuvqHvxDpxQM2mlY8iffxn58y8lvWx8ssMjtHUF2/94BeGtK8g//grKvvxr/NkDs5MLXRESERERkT63emczp/1hMVlpfl69RkVQV/G2Rupe+DW1L/wKp72JvGMuoOiM75I5/rh+1aV1xqjpjL3hXWqe+hm7n76F1lUvM+zKu8mZtiDZofUJXRESERERkYO2rqaFE3/7FgAL/2E+R5bmJDmi/iPe3kzdy7+hWzUu8QAAIABJREFU9vn/wGmtJ3fW+ZSedyMZo6YlO7T9at+0lB13XUF4x2oKTrqKoZf8En9mXrLD6jVdERIRERGRPrOpto1Tf7eYuLW8/i0VQR2ccBt1r/yW2ud+Qbx5NzkzzqH0/JvIrDg62aH1WuaY2Yy5aRk1f7uR2udup33DO1T85M0BVQztjwohERERETlg2+rbOfX3b9EWjfPat+YzuSw32SElnRMJUf/a79n97K3EG3eRPXUBQ86/icxxc5Id2kHxBTMYeuGtZE88ma2//hyVv72QUdc9g/EPjhJicGQhIiIiIodNVVOIU3+/mLq2KK9eM49pwwbPVYKDYWNRGt68l5onbyZWV0nW5FMZcu1jZB15XLJDOyRypi2g/PLfUXXPN9j5wLWUXfG7fvVs08FSISQiIiIivbazKcRnfr+Ync0hXvrmPGaNHJi9ih0K1nFoeucvVD9+PdHqDWQeMY/hV99P9qRTkh3aIVd48lVEqjdQ++ytBIeOp/jMf0p2SJ+aCiERERER6ZWd3pWgrfXtPP+NuRw7ujDZISWFtZaW95+i+rGfEK78kPRR0xl53dPkTD97UFwp6cmQL/2caPUGdj38A9JKx5A3+wvJDulTUSEkIiIiIvuVWAQ9d9VcThh7eL70sz+x1tK66mWqH/sJoY3vEiw7kuHf/gt5x1yA8fmSHV6fMz4fw66+j2jdNrb/4VLSCl8fsM8/AQz+T0xEREREPpWdTSFO+d2eIujEcalXBIW2LGfLL05j6+1nEG/cSfnX/4dxt6wif+5FKVEEdfAFMxn5/ScJ5Jex9defI1KzOdkhHbTU+dRERERE5IB1FEHbGlKzCIo2VLHjf77OxhuOJrxtJUO/cgfjbltL4YlXDpre0w5UIG8Io/7xOWwswrZfnU28tSHZIR0UFUIiIiIi0q3EIuj5b6RWEeRE2ql56ues/+F4GhY9QPGCf+KI29ZRfMZ38aWlJzu8pEsfNpGR332C8K51VN75JWwsmuyQDpgKIRERERH5hKouRVCqPBNkraVx8YOs/5cJ1Dz2E3KmfpYj/n0NQy++HX926vaQ153sSScz7Mo/0rr6Faru+xbW2mSHdEBS83qeiIiIiPSoqinEqSlYBLWtX8yuB6+jfcM7ZIw+muHffIDsiSclO6x+reD4K4hUb2D3kz8lOPQISs7512SH1GsqhERERESkUyoWQZHqjVQ/+iOa3nmYQEE5w666h/zjLk+pThA+jdLzbyJaswkTGFi3DKoQEhEREREA1u9u5aw/vsOOplBKFEGxllp2P/kz6l75LcafRsm5/4+Ss36ILyMn2aENKMYYhl19/4D7DiUVQiIiIiLCmxtrOe+eJQC8ePWxzB9TlOSI+o4Taafupd+w+5lbcNqbKTjx65SefyNphcOSHdqANdCKIFAhJCIiIpLy/ryskisfXkFFUSbPXjWXI0qykx1Sn7COQ+PiP1P96I+J1W0jZ8Y5DLngVjJGTEl2aJIEKoREREREUpS1lptfXMuNL67lpHHFPP7V2RRlBZMdVp9o+fAlqh/+IaGty8momMXwq+8je9IpyQ5LkkiFkIiIiEgKCsfiXPXICv532XaumD2Cuy6YTjAw+DoHCG9fzc6H/onWD/5OWkkFw7/1EHlzLlRHCKJCSERERCTV1LZGOP/eJbyxsY6fnTmBH31m/IB8xmNfYs27qXniRupf+z2+jByGXvwfFJ52rb4MVTqpEBIRERFJIWtrWjj77nfZ1tDOQ5cezcUzhyc7pEPKxiLUvXwnNU/ejBNqofCUayg9/0YCuSXJDk36mT4thIwxC4A7AD9wt7X21i7L/xG4CogBNcCV1totfRmTiIiISKp6ff1uvnDvUvw+w6vXzBtUPcNZa2l+70mqH/4BkV3ryZ66gLJLfkn68MnJDk36qT4rhIwxfuC3wOlAJbDEGPOUtXZ1QrP3gdnW2jZjzLeAXwAX9VVMIiIiIqnIWst//t9GfvjMGo4ozuLZq+YybhD1DBfaspydD/0jbWteIzhsEiP/8Tlyp5+Z7LBSynu7K0n3B5hSWJbsUHqtL68IzQHWW2s3Ahhj/gKcC3QWQtba1xLavw1c2ofxiIiIiKScplCUKx9ewWMrqzh/ahn3XDSD/My0ZId1SETrKql5/Hoa3rwXf3YRZZfdSeHJV2MCgyO/gSAUi3LT8pe4/cPXWTB8As+c/vVkh9RrfVkIDQe2JUxXAnP30f7rwPPdLTDGXA1cDTBq1KhDFZ+IiIjIoPZhVRNfuHcpG+vauP2cyfzTyWMHRacI8dYGdj97K3Uv3gHWoeiz11H6+Z/gzy5Mdmgp5a1dm7ly0cN83FjDlePn8MtjPpfskA5Iv+gswRhzKTAbOKm75dbau4C7AGbPnm0PY2giIiIiA9L/Lqvkm4+uJC89wKvXzOPEccXJDulTc6Jh6l/5Lbuf+jnxtnry532F0i/8lGBpRbJDSymt0TA/WvY8v1mziFE5Bbxwxjc4Y/iEZId1wPqyENoOjEyYHuHN24sx5jTgx8BJ1tpwH8YjIiIiMuiFY3Gue3IVv3trCyeOLeIvl82iPC8j2WF9KtZxaFz8IDWP/4To7i1kH3UGQy+8jYzRM5IdWsp5Zcc6vrHor2xqqePaScfx77POImeAdknel4XQEmC8MWYMbgF0MfDlxAbGmJnAH4AF1trqPoxFREREZNDbUtfGBfcvY8m2Bn5w8jhuOWsiAf/A/uLQlg9fovrhHxLaupyM0TMpv/JucqacluywUk5jpJ0fLHmGP659h/F5Jfzfmd/mhLKxyQ7rU+mzQshaGzPGXAu8gNt99p+stauMMTcDS621TwG3AznAX737Vbdaaz/fVzGJiIiIDFbPrt7F5Q+9T8yxPP7V2Zw/tTzZIX0qrR+/Qc3j19P20euklVQw/Jo/kzf3YoxvYBd2A421lqe3rebbix+nqr2JHxx1MjfN/CyZg6BDij59Rsha+xzwXJd51yeMq5wXERER+RTq2iJc9+Qq7l9aybTyPB69YhbjS3OSHdZBa1u/mJrHr6d11csE8ssY+pU7KDzlm/gG6O1XA9m6xhq+986TPL/9I44qKOOJU6/gmNLB03FZv+gsQUREREQO3BMfVPGtxz6gtjXC/zt9PD8+bTzpAX+ywzoo7RvepfqJG2j94O/4c0sZeskvKTzlGnzpWckOLeW0RMP8fMUr/GrVQtL9AX415/NcO+k40nwD89jqiQohERERkQGmujnMtU98wF9XVDFzeB5//8ZcZgzPT3ZYB6V983vUPHEDLcufwZ9TzJALb6PotH/Alz54vvB1oLDW8vCm5fzzkmfY3tbIFUfM5tZZZ1GWlZfs0PqECiERERGRAcJay0Pvb+e7T3xIczjOz8+cyA9OGUfaAOwQIbRlOTV/u4nm9/6GL7uQIV+6hcLTrsWfmZvs0FLSB3VVfOedJ1i4cyNHFw/nkZMvY/7QimSH1adUCImIiIgMANsb2/nWox/w9OpdzB1VwJ8umsHksoFXNIS2LKfmyZtpXvYEvqx8Ss+/iaIzvoc/a2Be0Rro6sNt3Pj+i/z2o7fID2bw+3lf5Koj5+JPgU4pVAiJiIiI9GNxx/LHt7fwr8+uIRxz+OXnJ/O9E8bi95lkh3ZAQltXuFeAOgqg8250C6DsgmSHlpLaYhHuXLOIW1e+SmM0xDcnHMtPZy6gOCN1bklUISQiIiLST72+fjff+9sqVlY1cfK4Yu66YNqA6xEutHWFewVo6eP4MvMoOe8Gis/4vgqgJIk5cf60bgk3LX+RHW1NnDViIrfMOovpRcOSHdphp0JIREREpJ/ZUtfGPz+9mkdXVjGqMJNHLp/Fl6aV433v4oAQ2rqSmidv6lIAfQ9/dmGyQ0tJjnV4bPMH/OS9v7O2qYZ5paN56KSvcGLZuGSHljQqhERERET6idZwjNteW8/tr23AGLh5wQT++eRxZKYNjG6L4+1NNL3zCA1v/In29YvdAujc6yn+7PdVACXRS9vX8m/LnmNZbSVTCoby5Ge+xudGTh5QhXVfUCEkIiIikmTWWh5evoMfPL2aysYQl8wczm1nT2JkYWayQ9sv6zi0fbyQhv+7h6alj2Ij7QTLJzLkwtsoPPkbKoCS6K1dm7n+/Rd4pWodo3MKue+Ei/nK2KNToiOE3lAhJCIiIpJEizfX8cNn1vDmpjpmDs/joUuP5vixxckOa78iu7fQ+OZ9NLx5L9GaTfgy88iffxkFJ3yNzHFzU/5qQ7JYa3lpx1puWfkKC3dupCQ9m1/POZdrJs4j3a9T/0R6N0RERESS4J0t9dzwwse88HENQ3KC/PGCaXxtzqh+3RucdRxaVjxL3Uu/oXX1y2At2ZM/w5Av/JTcWefjS89Kdogpy7EOf9uyiltWvsKy2kqGZ+Xz6znnctWRc8hOS092eP2SCiERERGRw2jptgZueOFjnltTTXFWGredPYl/OK6C7PT+e1oWb2+i4Y17qHvpN0SrNxAoGkHpuTeQf/wVBEsrkh1eSos6cR7a+D63rnyVNY3VHJFbwt3HXcCl42bpCtB+6N0REREROQzeq2zghhfW8szqXRRlpXHLWRO59rgx5Gb039OxyK711L30GxreuAcn1EzmEfMZ8qVbyJt1PiaQluzwUlpbLMK965bwiw9fZ0tLPdMKy/nLSZfypYppegaol/rvT56IiIjIILB8eyM3vvAxT67aRWFmGj87cwLfOX4MeRn9s5Cw1tK6+lXqXryDlhXPgC9A/pwLKTrje2SOPSbZ4aW8Tc21/PdHb/E/a9+lPtLOvNLR/PbY8zlrxCQ9l3WAVAiJiIiIHGLWWhZuqOX21zfw3JpqCjLTuHnBBL57/BjyM/tnARStq6Rx8Z9pePM+IjvW4M8tpeTzP6Hw1G+RVlCe7PBSmrWWl3es4841i3h622p8xvCF0UfxnUnHc/zQMSqADpIKIREREZFDJO5YHv+gittf28CSbQ2U5gS5eYF7BaigHxZATqiFpmVP0LjoflpXvwLWkjn+OIZddQ95cy/GF8xIdogprTka4v71y7hzzSI+aqymNCObH007lWsmzmNEdkGywxvwVAiJiIiIfEptkRj3LqnkVws3sKG2jfEl2fz+S1O5fPbIfvdlqNaJ07rmdRoX3U/T0sew4VbSSsdQcu71FMy/jODQcckOMeWtrNvB3Wvf5d71S2iOhjmmZCT3n3AJF1RMI0PPZh0yKoREREREDtLuljC/XbSZOxdtZndrhLmjCvjFOZM596iyftUNtrWW8LYPaHz7QRoX/5lYXSW+rHzyj/0yBcdfTub443R7VZLVhlp5aOP73LN+Ce/VbifN5+fCiul8Z/JxzC0dnezwBiUVQiIiIiIHwHEsr2+o5Z4lW3l0RRWhmMPnJg/lB6eM4/gxRf2qoIjUbKLx7YdoWvwg4e2rwOcnZ+pnyb/4l+TO/By+YGayQ0xpMSfOi9vXcs/6JTy1dRURJ87MouHcMfdcvjz2aEoyspMd4qCmQkhERESkF7bWt3HvkkruWbKVzXXt5GcE+OoxI/nO8WOYXJab7PA6xZqqaXrnERrffpD29YsByDzyeMou/2/yjvkSgbzSJEcoaxp2ce+6JTyw4T2q2psoSc/m2xPn89XxxzC9aFiyw0sZKoREREREehCKxvnbhzv507tbeXndbqyF08aXcMuZkzhvalm/ef4nWr+DlhXP0bT0UVpXvQxOnPRR0xly4a3kzb2YYIlurUq2tY01PLJpBX/dvIKV9VX4jY+zR0zia+OP4awREwnqy08PO73jIiIiIgkcx7J4Sz0PvredB9/fTkN7lNGFmdxw+pFcccxIKoqykh0i1onTvnEJLSuepWXFs4S2vA/gdnpw9r+Sd+wlZIyYkuQoZW1jDX/dvIJHNrnFD8BxQyr49ZxzuXjsDIZm9p8rialIhZCIiIikPGst729v5C/v7+DhFTvYWt9ORsDHF6eVc+WcUZw8rhhfkjs/iLc20PLhC7Qsf5aWD54n3rwbfH6yjpjPkAtvJWf62aQPn9KvnlFKRWsba3h080oe2byCFXU7gD3Fzxcrpqrb635EhZCIiIikrI92NfOX5Tt46P3trK1pJeAzfHZCKbecOZHPTykjNyN5p0o2FqV947u0rnqZllUv0b7hbXDi+LOLyJl2JjnTzyZn6mfx5xQlLUaBSDzGG7s28WzlGp7Ztpp1TbsBmK/ip99TISQiIiIpw1rLqp3NPL16F48s38HyHU0YAyePK+afTx7HF6eVU5QVTFpskaqPaV31Ei2rXqZtzWs4oWYwhoyK2ZSc/S/kTD+bzHFzMb7+8WxSqtrV3szzlR/xzLbVvLhjLc3RMOn+AKeUjeO7k47n3FFHMTJHxU9/p0JIREREBrVQNM7CDbU8s3oXz6zZxea6dgDmjCrgP8+dwoXThzEsP+Owx2WtJVqziba1b9D20UJaVr1ErK4SgLTSseTP+zLZU04je9KpuuqTZJF4jHdqtvJK1Tqer/yId3dvA2BYVh6XjJnJ2SMn8ZnyI8hOS09ypHIgVAiJiIjIoFPVFOK5NdU8s3oXL62toTUSJzPNx2njS/m3U8dz9uQhDM8/vN+hY5044W0fuIXP2jdpW/sGsQb3AXp/dhFZk08l5/Onkz3lNIJDxh7W2GRvccfh/brtvFq1nler1vPGro20xaIYDHNKR/LTmQs4Z+QkphcN0zNZA5gKIRERERnwalsjvLGxloUba1m4oZb3tzcBMLIgg8tnj+CcyUM55YiSw9rdtRNuo33TEtrWLaLt4zdoX/8WTrsbV6BoJFmTTiHryBPIOvJ40odNxvh8hy022ZtjHVbV72Lhzg28UrWe13duoCHiXjmcXDCUK8fP4TPl4zmpbCyF6cnvNVAODRVCIiIiMuBUN4f5P6/oWbixlg+qmgHICPiYV1HIz86cwOcmlzG1PPew/cc+Wr+D9vVv0bZ2EW3rFhHa+j7EYwCkD59C/rGXuIXPhBNIKx51WGKS7jVHQ7xbs41F1Zt4q3oLb9dsoTESAqAip5Avjp7KqeVHcGr5EZRl5SU5WukrKoRERESkX7PWsm53K4s31/PW5jre3FTH6l0tAGQF/cwfXciFC4Zx0rhi5owqID3Q91d9nFALocoPCG1+j/b1i2lbt4jo7s0AmLQMMsfOofjMH5A1fj6ZR8wjkFPc5zFJ96y1bGyu5Z2arbxVvZlF1ZtZWV+FYy0Gw1GFZVw8Zgbzh1Rw/NAxjM3VZ5UqVAiJiIhIv9ISjvHu1gYWb6lj8eZ63t5ST21bFID8jADzKgq5bNYIThpXzKwRBQQDfXdLmbWWWP0OQluXE9q6nPDWFYS2LidSvR6sBSBQUE7m+OMoOv27ZI0/jozRMzCB5PQ8l+qstWxormXp7m0sq61k2e5K3qvb3nm1JzsQ5NjS0fxk+mnMH1LB3JJRFKQf3mfFpP9QISQiIiJJ09geZcWOps7XssoGVlY14bg1BpOG5nDuUWXMG13IvIoiJg3J6bMvNnUi7YS3rya8bSWhbSsJV35AaOty4i21nW3SSseSMXoG+fMvJWPUDDJGzyBQNFIPzCdBOB7jo8ZqPqivYkXdjk8UPUGfn6mF5Vw8Zgazikcwu2QkUwvLCKjrcfGoEBIREZE+5ziWTXVtCUVPIyuqmjq7sgYoyQ4yc3gePz5tPPMripg7qoDCPvhOH+vEidZsJlT5gVvseEVPZOc6sA4AJphJ+vAp5B59XmfBkz5yGv5MPS9yuDnWYWNzHR/UV/Fh/c7O4dqm3cS9zyvo8zOtaE/RM6t4BEcVlhH061RXeqajQ0RERA6ZcCzOuppWPqpuYU11C2t2tbBmVzMf17TQHnVPWn0GJgzJ4dhRhXzz2NFMH5bH9GH5lOelH9IrK/HWesJVHxPZ+bE7rPqISNXHRKrXY2ORznZpQ8aRMXIaeXMuImPkNNJHTiU4ZJy+tPQwa46G+Lixho8ba1jbVOONV7O2qYa2WLSz3djcYqYWlvGF0VOZWljOUYVlHJlfSpo+LzlAKoRERETkgERiDlvq29hY28aG2jY21Layfncra3a1sLGujbh3X5sxMLowk0lDcjl1fAmTh+YyrTyPo8pzD0k31k40TLR2K9HdmxNeW4js3kxk1zriTdV7GvsDBIeMI1g2gZwZ5xAsn0D68ClkDJ+CLyPnU8civVMfbmNjcx2bWmrZ2FzHhqbazqKnyutaHMBnDBU5RUzIK+XksnEcVVjG1MJyJhcMJUdfWiqHiAohERER2YvjWHa1hNlS386WujY217ezsba1s+jZWt/e+QwPQGaaj7HF2UwflsfFM4cxcUgOk4bkMmFINlnBgz/ViLc1uoVO7VaitVu8omeLO757c+eXkXby+UkrGklaSQW5M84hWD6R9LIJBMsnECwdiwmkHXQssn/WWhoi7WxrbWBbawNbWxvY1FzHxuZaNrXUsbG5rvO7eToUp2dxZF4pnx1+JBPyhzAhv5QJ+aWMyy0hXbe1SR/TESYiIpJiWsIxtjeGvFe7W/DUt7Olvo3Nde1sbWgnHHP2WqckO8i44izmjy7isllZjCvOZmyxOzzQW9qcaJhY405ijTuJN+4k2lC1Z7x+u1vs1G3FaWvcaz0TCBIoGkla8Shypi4graQi4TWatMLhGJ089wnHOuwOtVLV3syOtkaq2po7C57Ewqc14ZZDcJ/dGZNbxNicYuaVjnbHc4sZm1vMmJwi8oIZScpIRIWQiIjIoGCtpaE9yq7mMNUtEXa1hKluDrOjKbSn6PHGm0KxT6w/NDed0YWZzByez3lHlTG6MJPRRVlUFGYyujCL3IyeTxms4xBvrSfWVE2suYZ4UzWxpmrizTXu0Jsfa9pFrHEnTmt9t9vx55YSKBxGWkkFWRNPIq14FGnFo91hyWgCeUMxvr7rKjvVdFzBqQ61UN3eQnWohZqQO6xqa6aqvYkdbU1UtTWxs72ZmHU+sY2yzFxGZhcwqWAoZwyfwKjsAkYmvMqzcvEZfWbSP6kQEhER6YciMYfatgi7WyPUtnpDb9qdF6W6JUx1S5hdzRFqWsNE4/YT2wn4DOV56QzPz2TSkBxOG1/K8PyMvV4jCzLJTPNj4zHibQ3EW+uIt1bjtNYT31xHdFU9NS217vzEYYs3bKvv/E6drnzZhQRyS/HnDSG9fCLZk04hkF9OIL+MQIE3zC8jkDdEt659CnHHoTEaoi7cRl24jdpwK7XhNmpDCePe/N2hVq/oaSXqxLvdXnF6FuWZeQzLymNywVCGZeVRnplHeVYuwzLz3WFWvm5fkwFNR6+IiMgh5jiW1kicpnCUxvYYTeEYje1RmsIxmkIxGkNRGtqj1LdFaQhFaWiP0dAe3fMKRWkJd3+CCpCbHqA4O43SbLfAmTksj/JsQ3l6jKFpEUrTohQFIhT4wuTaNmz7Tpz2RuJtDTitjcRrGr3pRpy2Birb6om31uMkPKzeHV9mHv7sIvw5xfizi0grGYM/Z890IG+Ie1Unbwj+vCEEckpU3PRCzInTEo3QFA3RFA3RGHGHTZFw57ymaJjGSIiGSDsNkXbqw+2d4w1e+54YDIXpmRSnZ1Gcns2I7HyOLh7BkIwchmTmuENvvDQjm5L0bHU7LSmhT49yY8wC4A7AD9xtrb21y/J04H5gFlALXGSt3dyXMYmIiFhriTmW9micUNQhFIvTHnVoi8Rpi8Y7h62R2F7zWiNxWiIxWsJxWsIxWiIdwz3zmsNu4dN5gcRa0oiRbiNkECFoI2R448XBOMVpcUrT4hwZiJEfiJGXESMvK0aOL0KeiZBNiEwbIsNpJ+iECMTaINKK09KKs7sFJ9SME26B+Cdvd2v2Xp2McYuZrAJ8mfn4M/MJFI0gfcRUt6DJKsSfU4TPG/qzC93CJ7sQf1ZhShY1jnUIxWO0x6LuMB6lPRZ1h3F3XlssQmssQlss6g33Hm/1Xi3RCC2xMC3RMC2xSOcw3M1n153ctHQKgpkUBjMpCGZSkVNEYbo7XhDMoDCY1VnwFHlFT3F6FgXBTPy6pVDkE/qsEDLG+IHfAqcDlcASY8xT1trVCc2+DtRba48wxlwM3AZc1FcxiYhI37HWEncscWuJxd1CIxp3vKE37TjE4olDSyTmEI07e8Ydh2jcXTfSOXQIR+NEYjFi0RjRWJRYLEY0GiUWixCNRt3xaIRYPEo8EiEWixGPhXG8tk4sQjwWxYlHsbEoPhsnQAy/dQgQI2DjpBFzXzbqDWMJwygZRCkycTJ9MTKJkWFipJsY6UQJ2ihBogRslKATxh+P4HfC+GJhDN3fNrZf/jR86dnuKyMHX9Abzy3FlFTgS89x52fk4svIwZ+R687LzN2zLDMff1Y+vqwCfOk5SX3GxlpL3DrEHIfYXsP4XtNRJ+69HGLWHUadODFnz7JIL4bheGyvYSQeI9w5jBGOxwk7MULxKKFYzBuPEY57QyfW6yKlK7/xkR0IkhVIIysQJDctnZxAkIJgJiOy8snxpt1hOjlpQfLSMsgLprvDLuM5aUE9ayNyiPXlFaE5wHpr7UYAY8xfgHOBxELoXOBGb/xR4E5jjLG2hxuN+5FwJMLDz9+X7DAOuQHw1h+kQ5BXL96bXu2lN+9xrz6HXsTj7N2m2zW6efj1wPcEtmM7nY3tXuuZhJzsJ0Y6JhO3YRMX7JlynM4mxttHx3hHI4vd+z30xt3N7h2ntbZzO53r2L2n3Z8L2xmXtXSe2Hbsy23q7AnUOgl52s5XR3uDdbfrrWc7cu5sB8bbnu1Yn09uB2u997YjDqdz/+4yd9vGOl4bx9uWu23TMd3R3mtjOuZhO8eNtz1Dx3KLj3jCfIvPOviw7ss6+HHn+62DD3e663w/Dn7rkI5DVsJyn7X4iXeOp9HzrWJ9xfH5cXxpWH8A63eHjj8NJxDE8QfdcX8ajj8LJ5BGzB8k7E+jORAk7g8STxx2GY/5g8Tg4JDgAAAO8klEQVS8YTzgjkc756cT86cR8wew1r0qYQHH++wd7xhJnO58YbFhByfciNPYgGO3dc532zh7tXW84sQddp3eM0xc1jHddXlHkRO3bkHjFsbeuNf2cAr6/KT7AwR9foK+AOn+jmGgc1mGP0BuIJ2MrEDntDtMI93nJzOQRqY/rXOY4U8jMxBw53nzswLBPUWP3x3q1jKR/q8vf0qHA9sSpiuBuT21sdbGjDGNQDGwuw/jOiSa2pqY/fjVyQ5DROQT4l6ZZo3BYnCMN43BGrwSxGCNwfHaObhdHztd1nG8ZR1tOpa7y3zEjTfPGOIJ68ZxpyPevLjx4Rg/cQLuNo3PjdMYYsbXuY7TMW7cdTrmxbx9xY3Pe+0ZdzBEfR1t3GHieOe8hDZR4yNm/MR8e5ZHjY+oz0/U+L3t+bEH0CV0Lz8ciDsQDgH7fqbDGPBh8BmDMQZf5zwfxrhtfMZ7JYybvcbdKxOd7fZq69urjd+b5/fa+Y2PgC/QuX7i/I72HeMdywM+f+f8gPHh93nDjvY+H2nGT8CbnzhM8/k710kzfnfa5yMtYVnnPOMn6BU1aT7fJ4ZpPnf9A+nSW0RSz4D4d4Ux5mrgaoBRo0YlORpXXlYez57/+2SHIZ7e/a07BH8QD9FtCaY3sfQm3N4k3qVNdycGvYtnT5vE1qa7NmbP0r3aJu7bG7XeiVu3a5i9x03naNf57gni3sH43MmO9RJi69iS+3EaOlY2gDG+vbbj8/kS1ndPSE1CLBiDz2cw+PZs2+fF09HeO/F05xsv/j37xRh3vwlx2sQ2Cfv2dtJlfV9nu+4+y64feWf+uA9wdqendfael9i++8+w61rdH397L+uc7ohzr8MmYdvdzO+6TuI7krj9vcYT1knczl7TZu/3reMY6Ph8E/dhuqzfOezSvqPA6WjfMS0iIodHXxZC24GRCdMjvHndtak0xgSAfNxOE/Zirb0LuAtg9uzZ/eLerfRgkK+e981khyEiIiIiIgehL5+6WwKMN8aMMcYEgYuBp7q0eQq4whv/EvDqQHg+SEREREREBrY+uyLkPfNzLfAC7t0Xf7LWrjLG3AwstdY+BfwP8IAxZj1Qh1ssiYiIiIiI9Kk+fUbIWvsc8FyXedcnjIeAC/oyBhERERERka7UIb2IiIiIiKQcFUIiIiIiIpJyVAiJiIiIiEjKUSEkIiIiIiIpR4WQiIiIiIikHBVCIiIiIiKSclQIiYiIiIhIylEhJCIiIiIiKUeFkIiIiIiIpBwVQiIiIiIiknJUCImIiIiISMpRISQiIiIiIinHWGuTHcMBMcbUAFv206wE2H0YwulvlHfqSdXclXdqUd6pRXmnllTNG1I398OV92hrbem+Ggy4Qqg3jDFLrbWzkx3H4aa8U0+q5q68U4vyTi3KO7Wkat6Qurn3p7x1a5yIiIiIiKQcFUIiIiIiIpJyBmshdFeyA0gS5Z16UjV35Z1alHdqUd6pJVXzhtTNvd/kPSifERIREREREdmXwXpFSEREREREpEcqhEREREREJPVYa/v8BYwEXgNWA6uA73nzi4CXgHXesNCb/xVgJfAB8BYwPWFbC4CPgfXAv+5jn1d4210HXOHNywKeBT7y4rh1H+v/HNgGtHSZfyLwHhADvpRCeY/ycnnfi/GswZL3vtoBXwVqgOXe66qB8Jl78/8OrPDi+D3g72H9bvcD/I+3/krgUSBnkOX9J6Aa+LDL/G5jToG8H2bPcb4ZWD5Y8u4p3i5t/gmwQMkgyjsDeDeh3U0Jy6719r3PnPtb3gnLn+p6DHdZ3tPvtTfYc5zvAP42yPLu9ud7oB7nvcm7p3i9ZRd48xxg9r6O80GY+wzgbe9YXwrMGQh5A69763f8nA7pYf1Z3v7XA//Fnsd9bgS2J6zf4/mqtfawFULlwNHeeC6wFpgM/KLjTQL+FbjNG5+f8GafCbzjjfuBDcBYIIj7y31yN/srAjZ6w0JvvBD3hPcUr00Q9xfimT3EfKwXd9eCoAKYBtzP/guhwZT3XcC3vPHJwObBkve+2uEWQncOtGPdW5bnDQ3wGHBxN+v3uJ+O9b3xX7HvX2gDKm9v+YnA0XyyIOg25sGed5c2vwSuHyx59xRvwvKRwAu4X9a9rxPEgZa3wfsHBpAGvAMc603PxP17tnlfOfe3vL3lXwAepOeTw97u5zHg8sGS9/5+vhmAx3kvP+8ef76BScAE3JPr3hRCgyn3F9lzLnMW8PpAyPsAPqt3cc9ZDfB8Qq43Av+8v/U7t9PbhofyBTwJnI5b8ZUnfAgfd9O2ENjujc8DXkhY9m/Av3WzziXAHxKm/wBc0k27O4Bv7CfWlh7m38t+CqHBlLe3rX9JiOetwZh313YcYCHUH3PHPQF6Griom/X3ux/cXzK/6/j8B0PeCW0q+GQhtN+YB2PeXT7vbcD4wZZ313gTph8FptOLomCg5o37D5/3gLld5h9QzsnOG8gB3sQ9SevpGO7N77U8oJ6Ef/gM9LwT1qvors1APM4PJO+u8XaZ9zq9OLkeTLnjFr0XJeznwQGS934/Ky+Wj7rbHgdYCB32Z4SMMRW4/4l6Bxhqra3yFu0EhnazytdxKz2A4bh/oDtUevO62m87Y0wB8DnglQNK4CANgrxvBC41xlQCzwHf6c1KAy3vHtp90Riz0hjzqDFm5L7W77KtCpKcuzHmBdzbJJpx/wge6Pr3ePFOBH7TzfqfMEDy3pfexPwJgyDvDicAu6y163rTeKDl3SVejDHn4v4RX7Gv9faznX6btzHGb4xZ7rV7yVr7zr7y2p9+kPdPca9Ytu0jzN7s5zzgFWtt0z6202mA5N2jAXycH1DeXX++P41BkPv3gduNMduA/8AtSg50O8k6b7vHGLPcGPP/jDGmh/Ur97H+td5525+MMYXdrN/psBZCxpgc3EvR3+/6y8e6ZZzt0v4U3Df4Xw5xHAHgIeC/rLUbD+W2e9jfYMj7EuBea+0I3EusDxhj9nn8DLS8e2j3NFBhrZ2Ge3/sfb3cZ7/I3Vr7Wdz/nKQDpx7E+l8DhgFrgIv2136w5J2wnU/E3J1BlvcluD8H+zXQ8u4arzEmC/gRcP2B7G8g5W2tjVtrZwAjgDnGmKMOdn/JztsYMwMYZ6194hBsbsAc558274F6nB9o3vuK9yD2PRhy/xZwnbV2JHAd7nO/B7Md4LD+bvuKtXYq7j/lTgAuO8D1fweMw31Gqgq3mOzRYSuEjDFpuG/un621j3uzdxljyr3l5bj/sepoPw24GzjXWlvrzd6Oe49rhxHAdmPMXK9yXG6M+XxP7RKm7wLWWWt/7e3Ln7D+zYcqZ2/bgyXvrwOPAFhrF+M+hFsyyPLeq52Xa621NuxN3o37cN4+9bPcsdaGcC9zn2uMGZmw/jW9XD8O/AX44iDKe196jHmQ593xz4Av4HacsL+2AyrvHuIdB4wBVhhjNnvbfc8YUzZY8k5o14D7MPSCnnLbl36S9zxgtvdZvQkcaYx5/UB/rxljSoA5uJ3kDKa8ezJQj/Ne591DvAdlEOV+BdAx/VfcY76/5421tmPYjPt81BzzyfO27d463a2/y/sHkAP8cX95H9C9kgf7wr3n/H7g113m387eD2H9whsfhdsLxPwu7QO4D1SNYc9DWFO62V8RsAn3vsVCb7zIW/Yz74P29TL2g35GaDDljXvZ86ve+CTc3nbMYMm7p3Z498Z64+cDbw+Ezxz33uLyhG09DFzbzfrd7sfL44iEnP4D+I/BknfCdir45DNC3cY82PP25i8AFvbi98OAyruneLtpt5l9P0Q+0PIuBQq88UzcjmDOOZCc+1PevTmGe7Mf4BrgvoFynPc27962GSjH+QF83vv9+ab3D+APmtxx7+Q42Rv/DLCsv+ftrV/itUnDveX3mh5i7tpZwlne/MTztuuAv+zzM9/fQXEoXsDxuJfTVpLQnR1QjPssxjrgZfacvN6N+xBjR9ulCds6C7c3iw3Aj/exzyu9D2k98DVv3ggvjjXspztk3J4yKnG7XKwEbvTmH+NNtwK1wKoUyXsysMg7qJcDZwyWvPfVDvh33K4kV+D+N3XiADnWhwJLvDg+xH2+J9DD+p/YD+7V4kW4XVN+CPyZfTxUPEDzfgj3snkU91j/uje/25gHe97esnvp4Y/OQM67p3i7abeZfZ8gDrS8p7HnKw8+JKEnQOC73ucfw/3H1t39Pe8uyyvY98l+j/vBPSleMFCO8wPMu8ef74F2nPc2757i9Zad770PYWAXCQ/yp0DuxwPLcM9f3gFm9fe8gWwv5pW451530PNXQczG/b22AbgTOrvPfgD33GUlbtfj5T3FYK3tXElERERERCRlHPZe40RERERERJJNhZCIiIiIiKQcFUIiIiIiIpJyVAiJiIiIiEjKUSEkIiIiIiIpJ5DsAERERLoyxsRxu0BNw+3i+X7gP637JXkiIiKfmgohERHpj9qttTMAjDFDcL9hPA+4IalRiYjIoKFb40REpF+z1lYDVwPXGleFMeYNY8x73ms+gDHmfmPMeR3rGWP+bIw51xgzxRjzrjFmuTFmpTFmfLJyERGR/kNfqCoiIv2OMabFWpvTZV4DMAFoBhxrbcgrah6y1s42xpwEXGetPc8Yk4/7jefjgf8E3rbW/tkYE8T9pvL2w5uRiIj0N7o1TkREBpo04E5jzAwgDhwJYK1daIz5b2NMKfBF4DFrbcwYsxj4sTFmBPC4tXZd0iIXEZF+Q7fGiYhIv2eMGYtb9FQD1wG7gOnAbCCY0PR+4FLga8CfAKy1DwKfB9qB54wxpx6+yEVEpL/SFSEREenXvCs8vwfutNZa77a3SmutY4y5AvAnNL8XeBfYaa1d7a0/Fthorf0vY8woYBrw6mFNQkRE+h0VQiIi0h9lGmOWs6f77AeAX3nL/ht4zBhzOfB3oLVjJWvtLmPMGuBvCdu6ELjMGBMFdgK3HIb4RUSkn1NnCSIiMmgYY7Jwv3/oaGttY7LjERGR/kvPCImIyKBgjDkNWAP8RkWQiIjsj64IiYiIiIhIytEVIRERERERSTkqhEREREREJOWoEBIRERERkZSjQkhERERERFKOCiEREREREUk5/x+MbXpiWL0qSgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 14min 48s\n"
     ]
    }
   ],
   "source": [
    "#get best parameters\n",
    "isolation_rate = study.best_params['isolation_rate']\n",
    "confirmation_rate = study.best_params['confirmation_rate']\n",
    "\n",
    "#reset timespan\n",
    "t_max = 80\n",
    "tspan = np.linspace(0, t_max, t_max+1)\n",
    "\n",
    "#reset differential equation: Fudan CCDC model\n",
    "def Fudan_CCDC(v, t):\n",
    "    #I0 must be positive\n",
    "    I0 = np.maximum(v(t)[0] - v(t)[1] - v(t)[2], 0)\n",
    "    #I0 at time s must be positive\n",
    "    Integral2, _ = integrate.quad(lambda s: f2(t - s) * np.maximum(v(s)[0] - v(s)[1] - v(s)[2], 0), 0, t)\n",
    "    Integral4, _ = integrate.quad(lambda s: f4(t - s) * np.maximum(v(s)[0] - v(s)[1] - v(s)[2], 0), 0, t)\n",
    "\n",
    "    dI = growth_rate * I0\n",
    "    dJ = growth_rate * Integral4\n",
    "    dG = lt(t) * (Integral2 - Integral4)\n",
    "    \n",
    "    return np.array([dI, dJ, dG])\n",
    "\n",
    "solution = ddeint(Fudan_CCDC, states_before_init, tspan)\n",
    "df_solution = pd.DataFrame({'Days': pd.date_range('2020-2-14', periods = t_max + 1, freq='D'),\n",
    "                            'Cumulative Infected': solution[:, 0],\n",
    "                            'Observed Infected': solution[:, 1],\n",
    "                            'Isolated Infected': solution[:, 2],\n",
    "                            'Potential Infected': np.maximum(solution[:, 0] - solution[:, 1] - solution[:, 2], 0)})\n",
    "df_solution['PCR Test Posives A'] = observed_J['A']\n",
    "df_solution['PCR Test Posives B'] = observed_J['B']\n",
    "df_solution['PCR Test'] = observed_J['PCR']\n",
    "\n",
    "#display table\n",
    "pd.options.display.precision = 0\n",
    "display(df_solution.head(25))\n",
    "\n",
    "#plot optimal solutions\n",
    "rcParams['figure.figsize'] = 14, 7\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(df_solution['Days'], df_solution.iloc[:,[1,2,4]])\n",
    "ax.legend(['Cumulative Infected', 'Observed Infected', 'Potential Infected'])\n",
    "ax.set_title('COVID-19 in Japan')\n",
    "ax.set_xlabel('Days')\n",
    "ax.set_ylabel('Number of I/J/I0 patients')\n",
    "ax.xaxis.set_major_locator(mdates.DayLocator(bymonthday=None, interval=7, tz=None))\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%Y-%m-%d\"))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}