fitting_models.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from datetime import datetime,timedelta\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import matplotlib.dates as mdates\n",
    "import matplotlib.ticker as ticker\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()\n",
    "\n",
    "import sklearn.metrics as sklm\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "# set font size\n",
    "# plt.rcParams['font.size'] = 24\n",
    "# plt.rcParams.update({'figure.figsize': (10, 7), 'figure.dpi': 120, 'font.size':10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expo_model(x,p,N0=3,X0=1):\n",
    "    return N0*(1+p)**(x-X0)\n",
    "\n",
    "# def logistic_model(x,a,b,c):\n",
    "#     return c/(1+np.exp(-(x-b)/a))\n",
    "\n",
    "def logistic_model(x,a,b,c):\n",
    "    return c/(1 + a * np.exp(-b * x))\n",
    "\n",
    "def gompertz_model(x,a,b,c):\n",
    "    return c*np.exp(-b*np.exp(-x/a)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustrate Exponential growth example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def f_exponential(t,r,x0):\n",
    "#     return x0*(1+r)**t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r = [2,1,0.5]\n",
    "\n",
    "# temp = []\n",
    "\n",
    "# for t in range(0, 31):\n",
    "#     r1 = f_exponential(t, r[0], 0.01)\n",
    "#     r2 = f_exponential(t, r[1], 0.01)\n",
    "#     r3 = f_exponential(t, r[2], 0.01)\n",
    "#     temp.append((t, r1, r2, r3))\n",
    "\n",
    "# exp_df = pd.DataFrame(temp, columns=['day','r=2','r=1','r=0'])\n",
    "# del temp\n",
    "# # print(exp_df)\n",
    "# # exp_df.style.hide_index()\n",
    "\n",
    "# format_dict = {'r=2':'${0:,.2f}', 'r=1': '${0:,.2f}', 'r=0': '${0:,.2f}'}\n",
    "# exp_df.style.format(format_dict).hide_index()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# temp = []\n",
    "\n",
    "# for t in range(0, 25):\n",
    "#     res = f_exponential(t, 1, 2)\n",
    "#     temp.append((t, res))\n",
    "\n",
    "# exp_df = pd.DataFrame(temp, columns=['hour','N'])\n",
    "# del temp\n",
    "# # print(exp_df)\n",
    "# # exp_df.style.hide_index()\n",
    "\n",
    "# format_dict = {'N':'{0:,.0f}'}\n",
    "# exp_df.style.format(format_dict).hide_index()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustrate Exponential, Logistic, Gompertz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1000.0, 150000.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABA8AAALlCAYAAACrT88FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAABYlAAAWJQFJUiTwAACoW0lEQVR4nOzdd3yV5f3/8dcVEsLeguACcWsd4KoTd90LZ2sdrda9tbYOwNmqddZVv1X0Z+vCOuu2oCiOiqvuCaiAKHuvXL8/7pPkJCSQcZKTc87r+Xicx7nXdd+fpNTkvHONEGNEkiRJkiSpNkXZLkCSJEmSJLVshgeSJEmSJGm5DA8kSZIkSdJyGR5IkiRJkqTlMjyQJEmSJEnLZXggSZIkSZKWy/BAkiRJkiQtl+GBJEmSJElaLsMDSZIkSZK0XIYHkiRJkiRpuQwPJEmSJEnSchkeSJIkSZKk5SrOdgFqOUII3wCdgHFZLkWSJEmSlHl9gVkxxn71bWh4oHSd2rZt22399dfvlu1CJEnKF/+b8j8WL10MwEY9N6J1q9ZZrkiSVKg++eQT5s+f36C2hgdKN2799dfvNnbs2GzXIUlS3uh0VScWL0rCg1d//yqd23TOckWSpEI1cOBA3nnnnXENaeucB5IkSU1kSdkSZi+aDUAg0LG0Y5YrkiSpYQwPJEmSmsishbMqtju36UxR8FcvSVJu8ieYJElSE5k+f3rFdpc2XbJXiCRJjWR4IEmS1ERmLJhRsd21TdfsFSJJUiMZHkiSJDWR9PDAngeSpFxmeCBJktREDA8kSfnC8ECSJKmJTF/gnAeSpPxQnO0CJEmS8pVzHkjNr6ysjGnTpjF79mwWLlxIjDHbJUlNIoRAaWkpHTt2pFu3bhQVNW3fAMMDSZKkJuKwBal5lZWV8e233zJv3rxslyI1uRgjCxYsYMGCBcydO5fVVlutSQMEwwNJkqQmYnggNa9p06Yxb948iouLWXnllWnfvn2T/zVWypaysjLmzp3L5MmTmTdvHtOmTaNHjx5N9jz/nyRJktREnPNAal6zZ88GYOWVV6Zjx44GB8prRUVFdOzYkZVXXhmo/PffZM9r0rtLkiQVMHseSM1r4cKFALRv3z7LlUjNp/zfe/m//6ZieCBJktREqkyY2NYJE6WmVj45oj0OVEhCCABNPjmo/6+SJElqIvY8kCQ1tfLwoKkZHkiSJDURwwNJUr4wPJAkSWoi0+c7YaIkKT8YHkiSJDWBhUsWMn/JfACKi4ppX+IEbpIK0/DhwwkhMHz48Aa1Hzp0KCEERo0aldG6VD+GB5IkSU1g5sKZFdtd2nRptjGpkgTJOPh8/+9OCIFBgwZlu4yCUZztAiRJkvKR8x1IUuLAAw9k6623pnfv3g1qf+qpp3L44Yez+uqrZ7gy1YfhgSRJUhNwvgNJSnTu3JnOnTs3uH2PHj3o0aNHBitSQzhsQZIkqQmk9zzo2qZr9gqRpBVYuHAhf/rTn9h4441p164dnTp1Yvvtt+ehhx6q8foYIzfeeCMbbLABbdq0YZVVVuHUU09l5syZ9O3bl759+1a5vrY5Dz744AOOOOII+vbtS2lpKSuttBIDBgzgzDPPZPHixRXXVZ/zoPx+AC+//HLFEI0QAkOHDq3yjLfeeovDDjuMVVZZhdLSUnr37s3uu+9e69em2tnzQJIkqQk4bEFSLli0aBF77LEHL7/8Muuttx6nnHIK8+bNY8SIERx22GG89957XHnllVXanHLKKdx222306dOHE044gdatW/PEE0/w1ltvsXjxYkpKSlb43A8++ICtttqKEAL77bcf/fr1Y9asWXz55ZfceuutXH755bXeZ9NNN2XIkCEMGzaMNdZYg2OOOabiXPocCHfeeScnnXQSrVq1Yr/99mPttddmypQpvP3229x6660ceuihDfqeFSrDA0mSpCZgeCApF/zlL3/h5ZdfZs899+SJJ56guDj5iDhkyBC23HJLrrrqKvbZZx+22WYbAEaPHs1tt93GOuusw5tvvkmXLl0AuPLKK9l1112ZOHEia6yxxgqfe88997BgwQIee+wx9t9//yrnpk+fTrt27Wptu+mmm7LpppsybNgw+vbtu0xvA4CPP/6Yk08+mU6dOjF69Gg23HDDKue/++67Fdaoqhy2IEmS1ASmL3DOA6mlCSF3Xs3lrrvuIoTAddddVxEcAPTs2ZOLL74YgP/7v/+rOH7PPfcAcOGFF1YEBwCtW7fmqquuqvfz27Ztu8yxrl27UlTUuI+qt912G0uWLOHiiy9eJjgAWHXVVRt1/0JkeCBJktQE7HkgqaWbPXs2X375JX369GG99dZb5vzOO+8MwLvvvltxrHx7u+22W+b6rbfeukoAsTyHHXYYrVq14oADDuDXv/419957L1999VVDvowavfHGGwDsueeeGbtnoTM8kCRJagJOmCippZs5cyZArUsolh+fMWPGMm169eq1zPWtWrWie/fudXr2lltuyejRo9l5550ZMWIERx99NGuttRbrrbce999/f32+jBqV17zKKqs0+l5KGB5IkiQ1AXseSC1PjLnzag7lyydOnjy5xvOTJk2qch1Ap06dAPjhhx+WuX7p0qVMnTq1zs//+c9/zlNPPcX06dN57bXXuPjii/nhhx848sgjefHFF+t8n5qUD6n4/vvvG3UfVTI8kCRJagKGB5Jauo4dO9K/f3++//57vvjii2XOjxw5EoABAwZUHNtss80AePXVV5e5/o033mDJkiX1rqO0tJRtttmGSy+9lJtuugmAxx9/fIXtioqKWLp0aY3ntt56awCeeeaZetejmhkeSJIkNQEnTJSUC4477jhijJx33nlVPoj/9NNPXHbZZRXXlPv1r38NwBVXXFExhAGSJR//+Mc/1vm5o0ePrtK+XHmPhuWttlCue/fufPvttzWeO+mkkyguLuayyy7j448/Xua8qy3Un0s1SpIkNYEqcx60dc4DSdlxzDHH1Hru1ltv5dxzz+WZZ57h8ccfZ5NNNmGvvfZi3rx5PPzww0yZMoXzzz+/yuSIO+64IyeccAJ/+9vf2HDDDTn44IMpKSnhySefpHPnzvTp06dOKyX85S9/4fnnn2fQoEGsueaadOjQgY8++ohnnnmGrl27csIJJ6zwHrvssgsPPPAA++67LwMHDqS4uJgddtiBHXbYgQ022IBbb72VE088kc0224z999+ftddem6lTp/L222/TsWPHip4VqhvDA0mSpCbgsAVJLUH50oo1ueGGG2jXrh0vvPAC1113Hf/85z+5+eabKS4uZpNNNuGGG27giCOOWKbdbbfdxnrrrccdd9zB7bffTvfu3TnwwAO58sorWXXVVenfv/8K6zr55JPp2rUrb775Jq+99hpLlixh1VVX5eSTT+acc85hjTXWWOE9brzxRkIIvPTSSzz99NOUlZUxZMgQdthhBwCOP/54NtpoI6699lpGjRrFY489Ro8ePdh444357W9/u8L7q6oQm2s2DrV4IYSxAwYMGDB27NhslyJJUk6LMdLmijYsWroIgPkXzqdNcZssVyXlv08++QSA9ddfP8uVFKYvvviCddZZh8MPPzwjKyao7ur6b3/gwIG8884778QYB9b3Gc55IEmSlGHzl8yvCA5KW5UaHEjKK5MnT6asrKzKsXnz5nHmmWcCcOCBB2ahKjU1hy1IkiRlmEMWJOWzG264gfvvv59BgwbRu3dvJk+ezEsvvcR3333HnnvuySGHHJLtEtUEDA8kSZIyzMkSJeWz3Xbbjffff5/nn3+eadOmUVxczDrrrMPpp5/OmWeeSQgh2yWqCeTcsIUQwuAQws0hhNEhhFkhhBhCuK8e7f+eahNDCGst57qjQwhvhRDmhBBmhhBGhRD2Wc71bUMIw0IIn4UQFoQQpoQQHgoh1DroJISwagjhrhDCxBDCwhDCuBDCDSGEWn/LCCFsE0J4OoQwLYQwL4TwQQjhzBBCq7p+DyRJUtOy54GkfLbLLrvwzDPPMGnSJBYuXMjcuXN59913Oe+88ygpKcl2eWoiORceABcBpwKbAt/Xp2EIYV/gOGDOCq67FhgO9AbuBO4DfgY8GUI4tYbrS4EXgEuAWcCNwIvAgcDbIYStamjTHxgLHAu8BVwPfA2cAbweQuheQ5v9gVeAHYBHgVuA1qm2D6zo65ckSc1j+vzpFduGB5KkfJCL4cFZwDpAJ+CkujYKIaxEEgQ8SPKhvbbrtgHOAb4CNo4xnhVjPAUYCEwDrg0h9K3W7GxgW2AEsFWM8fcxxiOBwUA74K4QQvXv9a1AT+D0GOMBMcYLYow7kwQB6wJXVKurU6r+pcCgGONvYoznkYQorwODQwiH1/X7IUmSmo49DyRJ+SbnwoMY48gY4xex/mtM/i31fsoKrjsx9X5FjLHizwYxxnEkf+kvJektAEBIBvSUtzk/xliW1uZxYDSwAbBjWps1gd2B8numGwLMBY4KIbRPOz4YWAl4IMb4dtozFpD0xoB6hCmSJKnpVJnzoI1zHkiScl/OhQcNEUI4BjgAODHGOHUFl++cen+2hnPPVLsGoD+wOvB5jPGbOrYp334+PWwAiDHOBl4j6bGwdR3regWYB2yTGkIhSZKyyJ4HkqR8k/erLYQQ1iCZg+C+GONjK7i2PbAKMCfGOKmGS75Iva+Tdmzd1Pvntdy2oW12T7V5aUVtYoxLQgjfABsCawKf1HJfAEIItQ3bWG957SRJUt0YHkiS8k1e9zxIzTNwD8kEiafXoUnn1PvMWs6XH++SI20kSVIWTF/ghImSpPyS7z0PziKZa2Dv9PkLMqA+8y2UL3LaYtrEGAfWeIOkR8KAejxTkiTVwDkPJEn5Jm97HoQQ1iZZseDuGOPTdWxW/tf7zrWcr+mv/ytq0ymLbSRJUhY4bEGSlG/yNjwgGf9fChwbQojpLypXPvgidewAgBjjXOB7oEMIoXcN91w79Z4+78Bnqfd1qFmTtwkhFAP9gCXA17XcU5IkNRPDA0lSvsnn8GAc8PdaXpNT1zyc2h+X1u4/qfdf1HDPPatdA/AVMAFYJ4TQr45tRqbed0/Ny1AhhNAR2BaYD7xRx7p2IFmdYUyMcWEN5yVJUjNyzgNJKlx9+/alb9++jb5PCIFBgwY1+j6ZkrfhQYzxvRjjb2t6UflX/D+mjr2X1vT21PuFIYSKQYohhL7AKcBC4O6058S0NlenhwEhhP2B7YGPgZfT2nwFPA+U3zPdMKA9cG+qJ0S5EcBPwOEhhM3TntEGuDy1e9vyvyuSJKk52PNAUraFEFb4GjVqVLbLzEmDBg0ihLDiC/NMzk2YmBpicEBqd+XU+89DCMNT2z/FGM9t6P1jjGNCCNcBZwMfhBBGAK2Bw4BuwGkxxnHVml0H7AMMBt4MIbwErA4cAswDjosxllVrczIwBrgphLALyfKKWwE7kQxXuLBaXbNCCMeThAijQggPANOA/UiWcRwBPNjQr1uSJGVGjNHwQFKLMWTIkFrPZeKv41rWSy+9lO0SmkTOhQfApsDR1Y6tmXoBjAcaHB4AxBjPCSF8AJwKnACUAe8A18QYn6rh+oUhhF2BC4AjSVZ5mAU8BgyJMX5cQ5uvUj0ILiUZirAXMAm4CRgWY5xWQ5vHQgg7kgQLBwNtgC9Jgo6bUr0gJElSFs1ZNIey1N8M2pe0p6RVSZYrklTIhg4dmu0SCk7//v2zXUKTyLlhCzHGoTHGsJxX3zrcY1Dq2i+Xc809McYtYoztY4wdY4w71hQcpF0/P8Y4JMa4doyxNMa4UozxkJqCg7Q238YYj40x9o4xto4xrhFjPKOm4CCtzWsxxr1ijF1jjG1jjD+LMV4fY1y6oq9bkiQ1Pec7kJRrvvnmG7p06UK3bt0YP358lXNz585l/fXXp1WrVrz8csVIbI455hhCCHz99ddcd911rLfeerRp04ZVV12Vs846i1mzZtX4rLFjx3LwwQfTs2dPSktLWWONNTj55JOZNGnSMteWP2PcuHHccccd/OxnP6NNmzb06tWLE044gZkza15o7rvvvuPUU09lzTXXpLS0lO7du7Pffvvx3//+d5lrhw4dWjGEY8SIEWy55Za0a9eObt26cfjhh/P9999XXDtu3DhCCBXfh/QhIOlzE9Q058HMmTO55ppr2HnnnVl11VVp3bo1K620Evvttx9vvPEGuSAXex5IkiS1WA5ZkJRr+vXrx//93/9xyCGHcMQRR/DKK69QXJx8VDz55JP59NNPGTp0KDvuuOMybc866yxeeeUVDj30UPbff3+ee+45brjhBkaPHs2rr75KmzZtKq596qmnOPjgg4kxMnjwYNZYYw3Gjh3LbbfdxuOPP85rr71W41CK888/n+eee459992X3XffnZEjR3LnnXfy5Zdf8p///KfKte+88w67774706ZNY4899uCggw7ip59+4rHHHmO77bbj0UcfZa+99lrmGbfeeitPPPEE++23HzvuuCNvvvkmDz74IO+//z7vvfcepaWldOnShSFDhjB8+HDGjx9fZUjIioaAfPLJJ1x44YXssMMO7L333nTt2pUJEybwxBNP8Mwzz/Dkk0/yi1/UNDd+y2F4IEmSlEHp4UHXtl1rv1CSmkFtwxbatGnDBRdcULE/ePBgTjrpJG677TYuvvhirrrqKu69917uvfdeBg0axMUXX1zjfV577TXee+891lhjDQCuuuoqDjnkEP71r39xzTXXVLSbM2cOxxxzDEuWLGHUqFFsv/32Fff485//zAUXXMAJJ5zA888/v8wz3njjDf73v/+x+uqrA7BkyRJ23nlnRo4cyVtvvcWWW25ZcfzQQw9lzpw5jBw5skrYMXHiRLbYYgt+85vfMG7cOEpLS6s849lnn+W///0vP/vZzyqOHXnkkdx///08/vjjHHrooXTp0oWhQ4cyatQoxo8fX68hIeuvvz4TJ06kR48eVY5/9913bLnllpx11lmGB5IkSYXEngdSyxWG5c4M+XFIZqYzGzZsWI3HO3fuXCU8ALjuuusYM2YMf/7zn1lllVW44IILWGmllfjHP/5BUVHNI97POOOMiuAAoKioiGuuuYbHHnuMu+66qyI8ePzxx5k6dSpHHHFEleAA4JxzzuH222/nhRdeYMKECRUhQblLLrmkyrHi4mKOPfZYRo8eXSU8+Pe//81XX33Fueeeu0wviT59+nD++edz5pln8tJLLy3T++D000+vEhwAHH/88dx///289dZbHHrooTV+/XXVuXPnGo+vuuqqDB48mJtvvrnGr70lMTyQJEnKoOnznfNAUstRnznV27Rpw4MPPsjmm2/OaaedRgiBESNG0KdPn1rb1DSUYc0112S11VZj3LhxzJgxgy5duvDOO+8AsPPOOy9zfXFxMTvssAPjxo3j3XffXeYD9Oabb75Mm9VWWw2A6dMr/5v7+uuvA9TaK+CLL74AkiEE1cODuj6jMV577TVuvPFGXn/9daZMmcKiRYuqnP/+++8NDyRJkgpFlZ4HpV2yVockNcQ666zDxhtvzJgxY9hggw3Yfffdl3t9r169ajy+8sorM378eGbOnEmXLl0qJjfs3bt3jdeXH58xY8Yy57p06bLMsfI5GZYurZw3furUqQA8/PDDy615zpw5DX5GQz366KMMHjyYNm3asNtuu9G/f3/at29PUVERo0aN4uWXX2bhwoWNfk5TMjyQJEnKIIctSC1XpoYC5LM//elPjBkzhh49evDRRx9x1VVXceGFF9Z6/Q8//MC66667zPHJkycDld31y9/Lj1dXvtpCbd3766K87eOPP85+++3X4Ps0hYsvvpjWrVvz9ttvs/7661c597vf/a7KShYtVc4t1ShJktSSOWGipFw1ZswYLrnkEtZdd10+/PBD1l13XYYMGcKrr75aa5uaPvR+/fXXfPvtt/Tt27fiL/qbbbYZAKNGjVrm+iVLllQ8Y8CAAQ2uf+uttwZg9OjRDb5HXbRq1QqoX4+EL7/8kg022GCZ4KCsrGy539+WxPBAkiQpg2YsnFGxbc8DSbli+vTpHHHEEbRq1YoHHniAXr168eCDD1JcXMwRRxxRMSSguhtvvJHx48dX7JeVlXHeeedRVlbGscceW3H8gAMOoFu3btx///288cYbVe5xww038PXXX7Prrrs2asz//vvvT//+/bnlllt4+umna7zm9ddfZ968eQ1+BkD37t0BmDBhQp3b9O3bly+++IKJEydWHIsxMmzYMD7++ONG1dNcHLYgSZKUQU6YKKklWd5yggcccACbbropAMcddxwTJkzgpptuqji2ySab8Je//IVTTz2VY489lieeeGKZe2y77bZsuummHHbYYXTu3JnnnnuO999/n4EDB3L++edXXNehQwfuuusuDjnkEHbccUcOOeQQVl99dcaOHcvzzz/PyiuvzB133NGor7WkpIR//etf7LHHHuy9995ss802bLrpprRr145vv/2W//73v3z99ddMmjSJdu3aNfg5u+yyCw8//DAHHXQQe+21F23btmWNNdbgqKOOqrXNWWedxYknnshmm23GwQcfTElJCa+99hoff/wx++67L08++WSD62kuhgeSJEkZ5JwHklqS2pZqhOSv4Ztuuik333wzjz32GPvttx+nnXZalWtOOeUUXnrpJR599FGuv/56zjrrrCrnr7/+eh599FHuvPNOxo0bR/fu3TnjjDO49NJLadOmTZVr999/f1577TWuvPJKnnvuOWbOnMnKK6/MiSeeyMUXX7zcVR3qauONN+b999/nuuuu46mnnuLuu++mqKiI3r17s9lmmzFs2DB69OjRqGf89re/Zfz48TzwwANcffXVLFmyhB133HG54cHvfvc7SktLueGGG7jnnnto27Yt22+/PXfffTePPPJIToQHoT5Ldyi/hRDGDhgwYMDYsWOzXYokSTlr09s35f0f3gfgnRPeYbPem2W5IqlwfPLJJwDLjCtX5h1zzDHcc889fPPNN/Tt2zfb5RS8uv7bHzhwIO+88847McaB9X2Gcx5IkiRlkD0PJEn5yPBAkiQpg6YvcM4DSVL+MTyQJEnKkKVlS5m1cFbFfqfSTlmsRpKkzDE8kCRJypD04KBzaWdaFbXKYjWS1HSGDx9OjNH5DgqI4YEkSVKGON+BJClfGR5IkiRliPMdSJLyleGBJElShtjzQJKUrwwPJEmSMsTwQJKUrwwPJEmSMiQ9POjatmv2CpEkKcMMDyRJkjKkSs+D0i5Zq0OSpEwzPJAkScqQ6fOdMFGSlJ8MDyRJkjLEOQ8kSfnK8ECSJClDZiycUbHtnAeSpHxieCBJkpQh9jyQJOUrwwNJkqQMcc4DSS3N559/ztlnn82AAQPo1q0bJSUldOvWja222opzzz2XsWPHZrvEFmfUqFGEEBg6dGi2S2lRDA8kSZIyxJ4HklqKGCPDhg1j/fXX5/rrryeEwGGHHcb555/Pr371K9q2bcvNN9/M5ptvzi233JLtcpUDirNdgCRJUr4wPJDUUlx66aUMHTqU1VZbjfvvv59tt912mWumTJnCDTfcwMyZM7NQoXKNPQ8kSZIyJD086NrGCRMlZcfXX3/N5ZdfTuvWrXnmmWdqDA4AevbsyZVXXsn5559f5fikSZM45ZRT6Nu3L61bt2allVbioIMOqnGIw/DhwwkhMHz4cF544QW23357OnTowEorrcSxxx7LjBkzAHj33XfZZ5996Nq1Kx06dGC//fZj3Lhxy9xv0KBBhBBYuHAhF110Ef369aO0tJT+/fszbNgwFi1aVOPX8umnn3LMMcew2mqrUVpaSq9evTjyyCP57LPPlrn2mGOOIYTA119/zc0338zGG29M27ZtGTRoEMcccww77bQTAMOGDSOEUPEaNWoU48aNq3Ksptfw4cOX879O7rLngSRJUgYsXrqYuYvnAlAUiujQukOWK5JUqO6++26WLFnCkUceyYYbbrjC64uLKz8WfvPNN2y33XZMnDiRnXfemSOOOIJvv/2Whx9+mH//+9888sgj7LPPPsvc44knnuCpp55in3324cQTT2TMmDEMHz6cb775hj/96U/ssssubL/99vzmN7/hf//7H08++SRfffUV//vf/ygqWvZv2oceeij//e9/GTx4MCUlJTz++OMMHTqUt99+myeeeIIQQsW1zz77LAcddBCLFy9m3333Za211uK7777jX//6F//+978ZOXIkAwYMWOYZZ5xxBqNHj2bvvfdmr732olWrVmyxxRYA3HPPPey4444MGjSo4vq+ffvSpUsXhgwZUuP38a9//StTp06lXbt2K/ye5yLDA0mSpAyoPmQh/RdbSWpOr732GgA777xzvdueeOKJTJw4kcsvv5wLL7yw4vjJJ5/MDjvswNFHH8348ePp0KFqQPrEE0/w0ksvseOOOwJQVlbGHnvswYsvvshee+3F3/72N375y19WXP+b3/yGu+66iyeffJL9999/mTo++eQTPvroI7p2TXpxXXHFFey000489dRT3HfffRx11FEATJ8+nSOOOIJ27drxyiuvsMEGG1Tc46OPPmKrrbbit7/9Le+8884yz3jnnXd499136devX5XjXbp04Z577mHQoEE1TppY07EhQ4YwdepUDjroIAYPHrzM+XzgsAVJkqQMcL4DKQeEkDuvRpg8eTIAq6yyyjLnxo0bx9ChQ6u8brjhBgC+++47nn/+eVZfffVlhjJss802HHHEEUybNo1//etfy9z3iCOOqAgOAIqKiio+4G+00UZVggOAX//61wC89957NX4NF198cUVwANCmTRuuuuoqAO66666K4/feey8zZsxg2LBhVYIDgA033JDjjz+ed999l48//niZZ5x//vnLBAcNce+993LppZey5ZZbct9999XYkyIf2PNAkiQpA5zvQFJLEWMEqLEH1Lhx4xg2bFiVY2ussQZnnnkm7777LgDbb789JSUly7Tdeeedue+++3j33XcrPvyX23zzzZe5vk+fPgAMHDhwmXPlwcZ3331X49eQHkSU23777SkuLq6oE+D1118H4P3336+xR8Dnn38OJD0ZqocLW265ZY3Pro+RI0fy29/+ln79+vHkk0/Stm3bRt+zpTI8kCRJygB7HkhqKXr37s2nn37K999/v8y5QYMGVYQLS5YsqRISlK+60Lt371rvC1RMgpiuc+fOyxwrn0theecWL15c47N69eq1zLFWrVrRvXt3pkyZUnFs6tSpANx555013qfcnDlzljm28sorL7fNinzyySccdNBBtG/fnn//+9/07NmzUfdr6fKzP4UkSVIzm75gesW24YHUQsWYO69GKF9d4aWXXqpXu/IP+eXDHqqbNGlSleua0g8//LDMsaVLlzJ16lQ6depUcay8lvfff58YY62vo48+epn7NWZumilTprD33nszd+5c/vWvf7H++us3+F65wvBAkiQpA+x5IKmlOOaYYyguLmbEiBF88skndW632WabAfDqq6+yZMmSZc6PHDkSoMaVCzLt5ZdfXubY6NGjWbJkSUWdAFtvvXXFuUxp1aoVkIQVNZk/fz777rsv33zzDf/3f/9XsbRjvjM8kCRJygDnPJDUUvTv35+LLrqIRYsWseeeezJmzJgar6s+/GDVVVdlt912Y9y4cRWTKJZ78803+ec//0nXrl058MADm6jySpdddhnTp1f26FqwYAF/+MMfADj22GMrjh977LF06dKFYcOG8dZbby1zn7KyMkaNGlWvZ3fv3h2ACRMm1Hi/X/3qV7z11lsMGTJkmbkf8plzHkiSJGWAPQ8ktSSXXHIJMUYuu+wytt12WwYOHMiWW25Jt27dmDFjBuPGjePFF18EYIcddqhod/vtt7Ptttty3nnn8fzzz7P55pvz7bff8vDDD1NUVMTdd99Nx44dm7z+9ddfnw033JDBgwdTUlLC448/zldffcXee+9dsYoDJB/0R4wYwYEHHsjWW2/NLrvswoYbbkhRURETJkzg9ddfZ+rUqSxYsKDOz1533XVZZZVVeOCBB2jdujWrr746IQSOOuoo3nzzTf71r39VBAw1TdJ4wAEHsOmmmzb2W9DiGB5IkiRlwPT5znkgqeUIITB06FCOOOIIbr/9dkaOHMk///lP5s6dS8eOHenfvz8nnXQSRx11VJVhCGuuuSZvv/02l19+OU8//TSjRo2iU6dO/OIXv+DCCy9kiy22aJb6H3roIS677DL+8Y9/MHHiRFZZZRWGDh3KBRdcsMxcBbvssgsffPAB1157Lc899xyjR4+mdevW9OnTh5133pmDDz64Xs9u1aoVjz76KBdccAEPPfQQs2fPJsbIdtttx7x584Bkosbqq1aU69u3b16GByE2cjIO5Y8QwtgBAwYMGDt2bLZLkSQp5xzxyBE88OEDANx34H38cuNfrqCFpEwrH99fCJPX5atBgwbx8ssv4+fU+qnrv/2BAwfyzjvvvBNjXHb9zBVwzgNJkqQMsOeBJCmfGR5IkiRlwE/zfqrYXqn9SlmsRJKkzDM8kCRJyoD08KBHux5ZrESSpMxzwkRJkqQMMDyQpMar77KKaj72PJAkSWqk+YvnM3fxXABKikro2LrplzGTJKk5GR5IkiQ10tT5Uyu2e7TrscwyYpIk5TrDA0mSpEZyyIIkKVuaa1lLwwNJkqRGMjyQWobyXj9lZWVZrkRqPuXhQVP3ejM8kCRJaiTDA6llKC0tBWDu3LlZrkRqPuX/3sv//TcVwwNJkqRGMjyQWoaOHZPJSidPnszs2bMpKytrti7dUnOKMVJWVsbs2bOZPHkyUPnvv6m4VKMkSVIjGR5ILUO3bt2YO3cu8+bN47vvvst2OVKzadeuHd26dWvSZxgeSJIkNZLhgdQyFBUVsdpqqzFt2jRmz57NwoUL7XmgvBVCoLS0lI4dO9KtWzeKipp2YIHhgSRJUiMZHkgtR1FRET169KBHD/+/KGWScx5IkiQ1kuGBJCnfGR5IkiQ1kuGBJCnfGR5IkiQ1kuGBJCnfGR5IkiQ1QozR8ECSlPcMDyRJkhph7uK5LFy6EIC2xW1pV9IuyxVJkpR5hgeSJEmNYK8DSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqBMMDSVIhMDyQJElqoBhjlfCge9vuWaxGkqSmY3ggSZLUQDMXzmRpXApAx9YdKS0uzXJFkiQ1DcMDSZKkBnLIgiSpUORceBBCGBxCuDmEMDqEMCuEEEMI99Vy7dohhN+HEP4TQvg2hLAohPBDCOHxEMJOK3jO0SGEt0IIc0IIM0MIo0II+yzn+rYhhGEhhM9CCAtCCFNCCA+FENZfTptVQwh3hRAmhhAWhhDGhRBuCCF0XU6bbUIIT4cQpoUQ5oUQPgghnBlCaLW8r0eSJGWe4YEkqVDkXHgAXAScCmwKfL+Cay8D/gT0Ap4G/gK8BuwN/CeEcHpNjUII1wLDgd7AncB9wM+AJ0MIp9ZwfSnwAnAJMAu4EXgROBB4O4SwVQ1t+gNjgWOBt4Drga+BM4DXQwjLDJoMIewPvALsADwK3AK0TrV9YAXfC0mSlGGGB5KkQlGc7QIa4CzgO+BLYEdg5HKufRb4c4zx3fSDIYQdST7sXxNCeDjGOCnt3DbAOcBXwBYxxump49eQfNi/NoTwVIxxXNotzwa2BUYAh8UYy1JtHgQeA+4KIfys/HjKrUBP4PQY481pz78u9TVeAZyYdrwTSZCxFBgUY3w7dfxi4D/A4BDC4TFGQwRJkpqJ4YEkqVDkXM+DGOPIGOMXMcZYh2uHVw8OUsdfBkaR/NV+m2qnyz+wX1EeHKTajCP5S38pSW8BAEIIIa3N+ekBQYzxcWA0sAFJ0FHeZk1gd6D8numGAHOBo0II7dOODwZWAh4oDw5Sz1hA0hsD4KTqX6skSWo6hgeSpEKRc+FBBi1OvS+pdnzn1PuzNbR5pto1AP2B1YHPY4zf1LFN+fbz1XojEGOcTTK0oh2wdR3regWYB2yTGkIhSZKageGBJKlQ5OKwhUYLIawB7ELygfuVtOPtgVWAOelDGdJ8kXpfJ+3Yuqn3z2t5XEPb7J5q89KK2sQYl4QQvgE2BNYEPqnlvgCEEMbWcmq95bWTJElVGR5IkgpFwYUHqb/M/4Nk+MH56UMTgM6p95m1NC8/3iVH2kiSpCZkeCBJKhQFFR6kljP8fySTGz4IXNvAW61wvoX0x7a0NjHGgTXeIOmRMKAez5QkqaAZHkiSCkXBzHmQCg7uAw4BHgJ+VcOki+V/ve9MzWr66/+K2nTKYhtJktSEDA8kSYWiIMKDEEIxcD9wOPBP4MgYY/WJEokxzgW+BzqEEHrXcKu1U+/p8w58lnpfh5o1eZvU19ePZPLHr2u5pyRJyjDDA0lSocj78CCE0BoYQdLj4F7gqBjj0uU0+U/q/Rc1nNuz2jUAXwETgHVCCP3q2GZk6n33EEKV/w1CCB1JhlXMB96oY107kKzOMCbGuLCG85IkKcOWli1l2vxpFfvd2nbLYjWSJDWtvA4PUpMjPgrsD/wdOLb60og1uD31fmEIoWvavfoCpwALgbvLj6eGPpS3uTo9DAgh7A9sD3wMvJzW5ivgeaD8numGAe2Be1M9IcqNAH4CDg8hbJ72jDbA5and21bwtUmSpAyZvmA6MTXVUNc2XSkuKqippCRJBSbnfsqFEA4ADkjtrpx6/3kIYXhq+6cY47mp7duBvUg+dH8PXBJC+byCFUbFGEeV78QYx4QQrgPOBj4IIYwAWgOHAd2A02KM46rd4zpgH2Aw8GYI4SVgdZLeDvOA42oILU4GxgA3hRB2IVlecStgJ5LhChemXxxjnBVCOJ4kRBgVQngAmAbsR7KM4wiSSSAlSVIz+HHujxXbDlmQJOW7nAsPgE2Bo6sdWzP1AhgPlIcH5cMIegCXLOeeo9J3YoznhBA+AE4FTgDKgHeAa2KMT1VvHGNcGELYFbgAOBI4C5gFPAYMiTF+XEObr1I9CC4lGYqwFzAJuAkYFmOcVkObx0IIO5IECwcDbYAvSYKOm2qYAFKSJDUR5zuQJBWSnAsPYoxDgaF1vHZQI55zD3BPPa6fDwxJvera5lvg2HrW9RpJ0CBJkrLI8ECSVEjyes4DSZKkpmJ4IEkqJIYHkiRJDWB4IEkqJIYHkiRJDWB4IEkqJIYHkiRJDfDTfMMDSVLhMDyQJElqAHseSJIKieGBJElSAxgeSJIKieGBJElSAxgeSJIKieGBJElSAxgeSJIKieGBJElSPS1auohZC2cBUBSK6NKmS3YLkiSpiRkeSJIk1dPUeVMrtru37U5R8FcqSVJ+8yedJElSPTlkQZJUaAwPJEmS6snwQJJUaAwPJEmS6snwQJJUaAwPJEmS6snwQJJUaAwPJEmS6snwQJJUaAwPJEmS6snwQJJUaAwPJEmS6umn+YYHkqTCYnggSZJUT/Y8kCQVGsMDSZKkejI8kCQVGsMDSZKkejI8kCQVGsMDSZKkejI8kCQVGsMDSZKkepi3eB7zFs8DoKSohI6tO2a5IkmSmp7hgSRJUj1MnTe1YrtHux6EELJYjSRJzcPwQJIkqR4csiBJKkSGB5IkSfVgeCBJKkSGB5IkSfVgeCBJKkSGB5IkSfVgeCBJKkSGB5IkSfVgeCBJKkSGB5IkSfVgeCBJKkSGB5IkSfXw03zDA0lS4TE8kCRJqgd7HkiSCpHhgSRJUj0YHkiSCpHhgSRJUj0YHkiSCpHhgSRJUh3FGA0PJEkFyfBAkiSpjuYsmsOipYsAaFvclnYl7bJckSRJzcPwQJIkqY7sdSBJKlSGB5IkSXVkeCBJKlSGB5IkSXVkeCBJKlSGB5IkSXVkeCBJKlSGB5IkSXVkeCBJKlSGB5IkSXVkeCBJKlSGB5IkSXVkeCBJKlSGB5IkSXX003zDA0lSYTI8kCRJqiN7HkiSCpXhgSRJUh0ZHkiSCpXhgSRJUh0ZHkiSCpXhgSRJUh2UxTKmzptasd+9bfcsViNJUvMyPJAkSaqDmQtmsjQuBaBj646UFpdmuSJJkpqP4YEkSVIdOGRBklTIDA8kSZLqwPBAklTIDA8kSZLqwPBAklTIDA8kSZLqwPBAklTIDA8kSZLqwPBAklTIDA8kSZLqwPBAklTIDA8kSZLqwPBAklTIDA8kSZLq4Kf5hgeSpMJleCBJklQH9jyQJBUywwNJkqQ6MDyQJBUywwNJkqQ6mDJ3SsW24YEkqdAYHkiSJK3A/MXzmbFgBgDFRcWGB5KkgmN4IEmStAKT50yu2O7VvhdFwV+hJEmFxZ98kiRJK5AeHqzcYeUsViJJUnYYHkiSJK3ApDmTKrZ7d+ydxUokScoOwwNJkqQVmDQ7LTzoYHggSSo8hgeSJEkr4LAFSVKhMzyQJElagSrDFux5IEkqQIYHkiRJK+CcB5KkQmd4IEmStAIOW5AkFTrDA0mSpBVwwkRJUqEzPJAkSVqOpWVL+WHuDxX79jyQJBUiwwNJkqTl+GneT5TFMgC6tulKaXFpliuSJKn55Vx4EEIYHEK4OYQwOoQwK4QQQwj3raDNNiGEp0MI00II80IIH4QQzgwhtFpOm6NDCG+FEOaEEGaGEEaFEPZZzvVtQwjDQgifhRAWhBCmhBAeCiGsv5w2q4YQ7gohTAwhLAwhjAsh3BBC6JrJr0WSJDWckyVKkpSD4QFwEXAqsCnw/YouDiHsD7wC7AA8CtwCtAauBx6opc21wHCgN3AncB/wM+DJEMKpNVxfCrwAXALMAm4EXgQOBN4OIWxVQ5v+wFjgWOCtVD1fA2cAr4cQumfia5EkSY3jfAeSJOVmeHAWsA7QCThpeReGEDqRfPhfCgyKMf4mxngeSfDwOjA4hHB4tTbbAOcAXwEbxxjPijGeAgwEpgHXhhD6VnvU2cC2wAhgqxjj72OMRwKDgXbAXSGE6t/rW4GewOkxxgNijBfEGHcmCQLWBa5o7NciSZIaz5UWJEnKwfAgxjgyxvhFjDHW4fLBwErAAzHGt9PusYCkBwMsG0CcmHq/IsY4Pa3NOJK/9JeS9BYAIIQQ0tqcH2NqUGTS5nFgNLABsGNamzWB3YHye6YbAswFjgohtG/k1yJJkhqpyrAFex5IkgpUcbYLaGI7p96freHcK8A8YJsQQmmMcWEd2jwDXJy6ZkjqWH9gdeDzGOM3tbTZPtVmZLVnPJ8eNgDEGGeHEF4jCRe2Bl5qxNdScL7+GqZMyXYVkqR88u6XleHBomm9eeONLBYjSWrxVloJ+vfPdhWZl+/hwbqp98+rn4gxLgkhfANsCKwJfJL6S/8qwJwY46TqbYAvUu/r1OUZjWyze6pNeXhQr6+llvsCEEIYW8up9ZbXrqW75ho4//xsVyFJyjuHTE5+wgI3XbkyN/0vu+VIklq2Y4+Fu+7KdhWZl3PDFuqpc+p9Zi3ny493aeD1Lb1NQblvuWtuSJLUQB3T/p4w22ELkqTClO89D1YkpN7rMn9Cuvpc35BnNGmbGOPAGm+Q9EgYUI9ntijz51dub7IJtGmTvVokSfnjvR6TKB8PuPGavWnbK6vlSJJauDXXzHYFTSPfw4Pyv8Z3ruV8p2rXrej6mv76X99nNGebgrJoUeX2Y49B377ZqkSSlC9ijHS4ajIsTvZf/vfKdDGcliQVoHwftvBZ6n2d6idCCMVAP2AJ8DVAjHEu8D3QIYRQU7/EtVPv6fMO1PqM5mpT09dSiNLDg9ats1eHJCl/zF40m3mL5wHQprgNnUtry/AlSQVr0SL46adk9vb334dx47JdUZPI954H/wF+CfwCuL/auR2AdsAr1VYn+A9wVKrN3dXa7Jl2TbmvgAnAOiGEfjWsuFBTm/JVF3YPIRSlr7gQQugIbAvMB9Lnc27I11JQ0sODkpLs1SFJyh+TZlddpjFZoVmSlBdihHnzYMaMqq+ZM2t/zZoFs2cn7+Wv9A8iAKeeCjff3PxfTxPL9/BgBPBn4PAQws0xxrcBQghtgMtT19xWrc3tJOHBhSGEx2KM01Nt+gKnAAtJCxVijDGEcDtwJXB1COGw8jAghLA/yTKNHwMvp7X5KoTwPMmKCqcA6f+yhgHtgTtSPSEa87UUlMWLK7fteSBJyoTJcyZXbK/cYeUsViJJqtWiRTB1avKaPh2mTUteNW3PmJFslwcF6R8iMmX27MzfswXIufAghHAAcEBqt/yn+M9DCMNT2z/FGM8FiDHOCiEcT/LBe1QI4QFgGrAfydKHI4AH0+8fYxwTQrgOOBv4IIQwAmgNHAZ0A06LMY6rVtZ1wD7AYODNEMJLwOrAIcA84Lj03gUpJwNjgJtCCLuQLK+4FbATyXCFC6vVVe+vpdA4bEGSlGmT5qT1POjoSguS1OSWLElCgB9/TIYC/Phj5eunn5Jz5e/lrzlzsltzcTF06lT56tMnu/U0kZwLD4BNgaOrHVsz9QIYD5xbfiLG+FgIYUeSD+MHA22AL0nCgZtijMusThBjPCeE8AFwKnACUAa8A1wTY3yqhusXhhB2BS4AjgTOAmYBjwFDYowf19DmqxDC5sClJEMR9gImATcBw2KM02poU++vpZAYHkiSMq36sAVJUgPMnw+TJ8MPP1S+pkxZdnvKlKRXQHN+rGnTBrp0qXx17rz8V6dOle/lr9JSKIBhbTkXHsQYhwJD69nmNZIP5/Vpcw9wTz2unw8MSb3q2uZb4Nh61lXvr6UQLF0KZam+HUVF0KpVduuRJOUHhy1I0nLMng3ffw8TJybvkyfDpEnJq3x78uRkroCm1KoVdOsG3bsn7+Wvrl2X3e/atWpY4PrudZZz4YFUE3sdSJKaQpVhC/Y8kFQoYkyGBnz7beXru+8qQ4LywKCpxvZ36wYrrbTsq0ePJCDo3r3qdqdOyV8Q1aQMD5QXDA8kSU3BOQ8k5aWFC5NAYNw4GD++8pUeFCxYkLnnlZTAyisnr549oVevZd979UoCgu7dkzkE1OL4v4rygss0SpKagsMWJOWkpUuTEODrrytf48ZVhgWTJmVmXoHSUlhlleTVpw/07l35Wnnlyu1u3QpiToB8Z3igvOAyjZKkpuCEiZJarAUL4Kuv4Isv4Msvk+30oGDJksbdv1MnWG21qq/0oGCVVZL5AwwFCobhgfKCwxYkSZm2aOkips6fCkBRKKJn+55ZrkhSwVmyBL75Bj77DD7/PAkKyl/fftvw3gNFRcmH/zXWSF59+8Lqqyev8qCgU6eMfinKfYYHyguGB5KkTPthzg8V2yu1W4lWRS7lI6mJzJ6dBASfflr5+uSTJCRI72JbHyuvDGuuWfnq27fyteqqjvVVvRkeKC8YHkiSMs3JEiVl3Lx5SSjw4YfJ66OPkvdvv63/vYqKkl4Da6+dvNZaKwkJ+vdPAoL27TNevgqb4YHyguGBJCnTnO9AUoOVlSXzELz/fvL63/+SoODrr+s/1KBPH1h33eRVHhSsvTb065dMWCg1E8MD5QXDA0lSprnSgqQ6mTMHPvigMigoDwvmzq37PYqLk0BgvfWS1/rrJ+/rruvcA2oxDA+UFwwPJEmZVmXYgj0PJEESFLz7LowdW/n69NO69yYoKkpCgo02gg03rHxfe23nIFCLZ3igvJAeHvjfXUlSJlQZtuCcB1LhWbAgCQreegvefrv+QUGvXrDJJslr443hZz9LehK0adO0dUtNxPBAeSF9Elp7HkiSMmHyXIctSAUjxmQ+gjfegDffTN7fe69uKx0UFSVDDDbZBDbdtDIwWNn/bii/GB4oLzhsQZKUaU6YKOWxBQvgv/+FV1+F115LAoOfflpxu6Ii2GADGDiw8rXJJq5soIJgeKC8YHggSco0l2qU8sjUqTBmTBIWvPpqMgwh/RfI2qy9Nmy9NWyxhUGBCp7hgfKC4YEkKZPKYhk/zPmhYt9hC1KOmTIFRo2CkSPhlVfg449X3KZrV9hqq8rXlltC9+5NXqqUKwwPlBcMDyRJmTRt/jQWlyVjnTuVdqJdSbssVyRpuaZOhZdfTsKCkSPho49W3Ga99WC77WDbbWGbbZJeBiE0fa1SjjI8UF4wPJAkZZLzHUgt3Lx5MHo0PP88vPQSvP/+8q8vLk6GHWy3XWVgsNJKzVOrlCcMD5QXDA8kSZk0eY4rLUgtSozwv/8lYcFzzyXBwcKFtV9fUpIMPdhpJxg0KJm3oJ09iKTGMDxQXkgPD0pKsleHJCk/OFmi1AJMnZoEBc89l4QGkyfXfm2rVsmkhjvtlLy22caJDaUMMzxQXkhfgteeB5KkxnLYgpQFMSYTGz71VPIaMwbKymq/foMNYPfdYbfdYPvtoWPH5qtVKkCGB8oLDluQJGWSwxakZrJgQTLRYXlgMG5c7dd265YEBeWBwWqrNVuZkgwPlCcMDyRJmVRl2II9D6TMmj0b/v1v+Ne/4OmnYe7cmq8LIZmrYK+9YI89YMCAZHiCpKwwPFBeMDyQJGWScx5IGTZ1KjzxRBIYPP981V/e0nXsCL/4BeyzD+y5pysiSC2I4YHyguGBJCmTHLYgZcAPP8AjjySBwahRsHRpzdf17w/77pu8ttvOX+akFsrwQHnB8ECSlElOmCg10PTpSVhw//0wcmTtEx5uuikcdFDy2mCDZIiCpBbN8EB5wfBAkpQpcxfNZfai2QCUFJXQrW23LFcktXBz58KTTyaBwTPPVF0GK93Pf56EBQcemPQ2kJRTDA+UF9J/RpWUZK8OSVLuqz5kIfgXUWlZS5bAc8/BffclcxnMm7fsNSEkSygeckgSGKyySvPXKSljDA+UF+x5IEnKFCdLlJbjww/hnnvg//2/ZE6Dmmy+ORxxBBx6KKy6avPWJ6nJGB4oLxgeSJIyxfkOpGqmToUHHoDhw+Htt2u+Zv31k8Dg8MNh7bWbtTxJzcPwQHnB8ECSlCmutCCRrIzw3HNw993JsISallZceWX41a+S18YbO+mhlOcMD5QXDA8kSZlSZdiCPQ9UaCZNgr//He68EyZMWPZ869ZwwAFwzDGw225Q7McJqVD4/3blBcMDSVKmOOeBCk5ZGbz4ItxxBzz+eNLroLqttkoCg8MOg65dm71ESdlneKC8YHggScoUhy2oYEyZkgxL+Nvf4Ouvlz3fo0cSGBx3XDKngaSCZnigvJAeHrhUoySpMZwwUXnvnXfgxhuTSRBrmstgxx3hd7+Dgw6C0tLmr09Si2R4oLyweHHltj0PJEmN4bAF5aUlS5IhCTfeCKNHL3u+S5ekl8EJJ9jLQFKNDA+UFxy2IEnKhCVlS/hx7o8V+z3b98xiNVIGzJgB//d/8Ne/wvjxy57fais4+WQ45BBo27bZy5OUOwwPlBcMDyRJmTBl7hQiEYAe7XrQupU/VJSjvvoKrr8ehg+HuXOrnisuhkMPhTPOgC23zEp5knKP4YHyguGBJCkTnO9AOe+99+DPf4aHHkpWUUjXvTuceCKcdBKsskpWypOUuwwPlBcMDyRJmZC+0oLzHShnxAivvAJ/+hM8++yy5zfaCM48E4480qEJkhrM8EB5wfBAkpQJ6ZMlukyjWryyMnjqqSQ0eP31Zc/vthucfz7ssguE0Pz1ScorhgfKC4YHkqRMqNLzwGELaqmWLoUHH4QrroCPP656LgQYPBh+/3sYODA79UnKS4YHygvpSzWWlGSvDklSbnPOA7VoZWXJXAaXXgqffFL1XOvWcPTRcO65sM462alPUl4zPFDOi7FqzwPDA0lSQzlsQS1SWRk88ggMGwYffVT1XIcOySSIZ50Fffpkpz5JBcHwQDlvyZLK7eJiKCrKXi2SpNzmhIlqUcrK4LHHYOhQ+N//qp7r2DGZBPGss6Br1ywUJ6nQGB4o5znfgSQpU9J7HjhsQVkTYzIR4sUXw/vvVz3XoQOcfjqccw5065ad+iQVJMMD5TzDA0lSJsQYq8x54LAFZcUbbyQrJIweXfV4+/Zw6qnJnAY9emSnNkkFzfBAOc/wQJKUCTMXzmTh0oUAtC9pT8fSjlmuSAXl88/hj39M5jZI17ZtEhqcdx6stFJ2apMkDA+UBwwPJEmZUGWlBec7UHP54Ydk9YQ77kiWYCxXXJxMhHjRRdCrV/bqk6QUwwPlPMMDSVImuNKCmtWcOfCXv8C11ybb6Q49FK64AtZaKzu1SVINDA+U8xYvrtx2mUZJUkNVWWnByRLVVMrK4P/9P7jgApg8ueq5HXeEq6+GLbfMTm2StByGB8p59jyQJGVClWELhgdqCm+9Baedlryn23BD+POfYa+9IITs1CZJK1CU7QKkxjI8kCRlgsMW1GQmT4Zjj4WttqoaHPTpA3fdlSzHuPfeBgeSWjR7HijnGR5IkjKhyrAFJ0xUJixaBDfdlEyIOHt25fHWreGcc5LVFTp0yF59klQPhgfKeYYHkqRMSO954LAFNdozz8CZZyZLMKbbbz+47jro3z8rZUlSQxkeKOcZHkiSMiF9zgOHLajBJk6E00+HRx6peny99eDGG2H33bNTlyQ1knMeKOcZHkiSMsFhC2qUpUvhlluSkCA9OOjUCa6/Hj74wOBAUk6z54FyXnp44FKNkqSGWLBkAdMXTAegVWhFj3Y9slyRcsr778Pvfgdvvln1+DHHJKso9OyZlbIkKZPseaCct3hx5bY9DyRJDZE+ZKFXh14UBX9FUh3MnQvnnw8DB1YNDtZdF0aOhLvvNjiQlDfseaCc57AFSVJjjZ85vmJ7jc5rZLES5YxnnoGTT4Zx4yqPtW6drKBwwQVQWpq10iSpKRgeKOcZHkiSGmvCzAkV26t3Xj2LlajFmzYNzjgD7ruv6vEdd4Tbb0/mPJCkPGSfPOU8wwNJUmMZHqhOnnoKNtqoanDQrVsyPGHkSIMDSXnN8EA5z/BAktRY42c4bEHLMWMGHHss7LsvTKqcH4Nf/hI+/TSZGDGEbFUnSc3CYQvKeYYHkqTGmjDLngeqxbPPwm9/C99/X3msVy+44w7Yf//s1SVJzcyeB8p5hgeSpMaq0vOgiz0PBMyaBccfD3vuWTU4OOII+OgjgwNJBceeB8p56Us1lpRkrw5JUm6KMTrngar6z3+SoQjfflt5rEePZELEgw/OWlmSlE32PFDOs+eBJKkxps6fyvwl8wHo2LojnUs7Z7kiZc2iRckyi7vuWjU4OPjgpLeBwYGkAmbPA+U8wwNJUmNUH7IQnPiuMH3xBRx5JLz9duWxbt3g1lvh0EOdEFFSwbPngXKe4YEkqTEcslDgYoR77oHNNqsaHOy2G3z4IRx2mMGBJGF4oDxgeCBJaozxM12msWDNmJH0NjjmGJg7NzlWUgLXXpusstC7dzark6QWxWELynmGB5KkxrDnQYF67TX45S9hfGV4xDrrwP33w4AB2atLklooex4o5xkeSJIaw/CgwJSVwWWXwQ47VA0OfvtbeOcdgwNJqkXBhAchhL1DCM+HEL4LIcwPIXwdQng4hPDzWq7fJoTwdAhhWghhXgjhgxDCmSGEVst5xtEhhLdCCHNCCDNDCKNCCPss5/q2IYRhIYTPQggLQghTQggPhRDWX06bVUMId4UQJoYQFoYQxoUQbgghdK3fdyR/GB5IkhrDYQsF5KefYK+94JJLkhABoEsXePhhuPNOaN8+q+VJUktWEOFBCOHPwFPAAOBZ4EbgHWB/4LUQwq+qXb8/8AqwA/AocAvQGrgeeKCWZ1wLDAd6A3cC9wE/A54MIZxaw/WlwAvAJcCsVE0vAgcCb4cQtqqhTX9gLHAs8Faqnq+BM4DXQwjd6/gtySuLF1dul5Rkrw5JUm6y50GBePPNpFfBc89VHtthB/jgAxg8OHt1SVKOyPs5D0IIKwPnAj8AG8cYp6Sd2wn4D3ApyYd9QgidSD78LwUGxRjfTh2/OHXt4BDC4THGB9Lusw1wDvAVsEWMcXrq+DUkH/avDSE8FWMcl1ba2cC2wAjgsBhjWarNg8BjwF0hhJ+VH0+5FegJnB5jvDnt+dcBZwFXACc2/LuVm+x5IElqqPmL5zNlbvKrQavQij4d+2S5ImVcjHDLLXD22VX/4vDHP8Kll0KrWjuVSpLSFELPgzVIvs4304MDgBjjSGA2sFLa4cGp/QfKg4PUtQuAi1K7J1V7RvkH9ivKg4NUm3EkvRZKSXoLABCSBaTL25yfHhDEGB8HRgMbADumtVkT2B0ov2e6IcBc4KgQQsH1tzM8kCQ11Lezvq3YXrXTqrQq8oNkXpkzJ1lN4bTTKoODLl3gySfhiisMDiSpHpolPAghrBdCOCuE8LsQQufmeGaaL4BFwJYhhB7V6toB6EgyXKDczqn3Z2u41yvAPGCb1LCDurR5pto1AP2B1YHPY4zf1LFN+fbz1XojEGOcDbwGtAO2ruF+ec3wQJLUUA5ZyGMffwxbbAEPpI04HTAgmRRxn1qnpJIk1SKj4UEI4ZIQwqQQQre0Y7sC7wLXknS7f6c5x+bHGKcBvwd6AR+HEP4WQrgqhPAQ8DzJvAO/S2uybur98xrutQT4hmS4x5oAqb/0rwLMiTFOqqGEL1Lv69TlGRluUxAMDyRJDTV+RtpkiV2cLDFv3H8/bLklfPpp5bETTkiWZ+zXL3t1SVIOy/ScB3sCn6Y+sJe7CogkXetXBk4mmeDvkgw/u1YxxhtCCOOAu4Dj0059CQyvNpyhvGfEzFpuV368SwOvb842NQohjK3l1HoratsSGR5IkhqqSs+DTvY8yHlLl8If/gDXXFN5rG1buP12+PWvs1eXJOWBTA9b6At8Ur4TQlgFGAjcGmO8PMZ4Ksmkgwdk+LnLFUI4n2RiwuEkQwbap+r6GvhHCOHq+twu9R7rWUZ9rm/IMxpaV84zPJAkNdSEWZXhgT0PctzMmbDfflWDg7XXTlZZMDiQpEbLdM+DrkB6r4NtST7MPpV2bCxVhwk0qRDCIODPwKMxxrPTTr0TQjiQZBjAOSGE22OMX1P5F/za5mbolHqfWe29tutr6jFQ32c0tE2NYowDazqe6pEwYEXtWxrDA0lSQ6UPW3DOgxz2xRdJcJA+TGGffeC++6Bzc0+3JUn5KdM9D34kGf9fbidgMfBm2rHWTfDc5SmfEWdk9RMxxnnAW6l6Nksd/iz1vszcASGEYqAfsISk1wIxxrnA90CHEELvGp6/duo9fa6CWp+R4TYFIX3VpZKS7NUhSco9TpiYB154Ydn5Df7wB3jsMYMDScqgTH+Ifw/YL4SwUQhhLeAw4NUY4/y0a/oCNU0s2FTKV0VYqZbz5cfL/379n9T7L2q4dgeSFQ3GxBgXph1fXps9q10D8BUwAVgnhFDTrD01tSkPP3YPIVT53y2E0JGkl8d84I0a7pfX7HkgSWqIslhWZalGw4McEyPceCP84hcwY0ZyrE0b+Mc/4MorXYZRkjIs0+HB1STd6t8n+Ut5Z+Av5SdDCG2AQcDbGX7u8oxOvZ+QmoOhQghhT5IP3QuAManDI4CfgMNDCJunXdsGuDy1e1u1Z9yeer8whNA1rU1f4BRgIXB3+fEYY0xrc3V6GBBC2B/YHvgYeDmtzVckq0OU3zPdMJJ5HO5N9YQoKIYHkqSG+GHODyxamvwQ6da2Gx1ad8hyRaqzhQvh+OPhzDOhLLWCdZ8+8MorcOSRWS1NkvJVRuc8iDGODiHsQ7KiQQT+EWN8Ju2SbYBxwKOZfO4KjABeBHYFPgkhPApMBtYnGdIQgAtijFMBYoyzQgjHp9qNCiE8QDKPw34kyyWOAB5Mf0CMcUwI4TrgbOCDEMIIkuEZhwHdgNNijOOq1XVd6vmDgTdDCC8BqwOHAPOA42KMZdXanEwSctwUQtiFZHLKrUiGh3wOXNjQb1IuMzyQJDVE+pCFNTo7WWLOmDIFDjooWXax3FZbwaOPQu+aRpBKkjIh0xMmEmN8Fni2lnP/oXJugWYRYywLIexF8tf6w4EDSYYeTAOeBm6KMT5frc1jIYQdST6MHwy0IVnW8ezU9cusaBBjPCeE8AFwKnACUAa8A1wTY3yqhusXhhB2BS4AjgTOAmYBjwFDYowf19Dmq1RviEtJhkjsRTIE5CZgWLUlMguG4YEkqSHGz3SyxJzz+eew557w9deVx446Cv72t2TIgiSpyWQ8PGiJYoyLgRtSr7q2eY3kw3l9nnMPcE89rp8PDEm96trmW+DY+tSV7wwPJEkN4WSJOWbMmGRFhalTk/0Q4Oqr4Zxzkm1JUpPK+KoHIYSiEMJpIYQ3QggzQwhL0s5tFkK4NYRQ24oBUr2UlcHSpcl2CM6NJEmqO4ct5JBHHoGdd64MDtq2TYYpnHuuwYEkNZOMhgchhNbACyR/4e8PzCaZU6DcN8BxwC8z+VwVrurLNPr7gySprhy2kANihOuvh0MOSSZJBFhpJRg5EvbfP7u1SVKByXTPg/NIJu8bBvQC/i/9ZIxxBvAKsEeGn6sC5ZAFSVJDVel50MWeBy3O0qVwxhlw9tlJiACw9trw+uvJBImSpGaV6fDgl8BrMcZLUysFLDOxIEnvA+N9ZYThgSSpocbPsOdBizVvHgweDDffXHls222T4KB//+zVJUkFLNPhQT/gjRVcM41k+UKp0QwPJEkNMXvhbKYvmA5A61at6dm+Z5YrUoUff0zmN3jsscpjhxwCL74I3btnrSxJKnSZDg/mA11WcM3qwIwMP1cFyvBAktQQ3876tmJ79c6rUxQyPoe0GmLCBNhuO3jzzcpj554LDzzgUoySlGWZXqrxPWD3EELrGOOi6idDCJ1J5jsYk+HnqkAZHkiSGsIhCy3QZ5/BbrvBt6lgp6gIbrwRTj01u3VJkoDM9zy4E1gN+EcIoVP6iRBCF2A40BW4PcPPVYEyPJAkNYTLNLYw77yT9DgoDw5at4aHHjI4kKQWJKM9D2KM94cQdgWOBfYDpgOEEN4GNgRKgVtijE9n8rkqXIYHkqSGcJnGFuSVV2CffWD27GS/fXt49NGkF4IkqcXI+AC/GONvgOOAj4GVgAAMAL4EfhNjPC3Tz1ThWry4crukJHt1SJJyS3rPA8ODLPr3v2GPPSqDg65dk4kRDQ4kqcXJ9JwHAMQYhwPDQwhtSYYpzIwxzm2KZ6mw2fNAktQQDltoAf75Tzj6aFiyJNlfeWV4/nn42c+yW5ckqUZNEh6UizHOJ1mBQWoShgeSpIZw2EKW3XprMp9BjMl+v35Jj4M118xuXZKkWmV02EIIoWsIYYMQQmm148eGEB4PIfwzhLBlJp+pwmZ4IEmqryVlS/h+1vcV+6t1Xi2L1RSgq66CU06pDA422ghefdXgQJJauEz3PLgS+BXQs/xACOE04AaSuQ8ADgghbB5j/DjDz1YBMjyQJNXXxNkTWRqXAtCrfS/aFLfJckUFZNgwGDq0cn+rreDpp6Fbt6yVJEmqm0xPmLgt8FJquEK5c4HvgR2AQ1PHzs7wc1WgDA8kSfXlZIlZECMMGVI1ONhll2SogsGBJOWETPc8WAV4qXwnhLABsBrw+xjjq6ljh5AECVKjGR5IkuqrymSJXZwsscnFCBdfDFdcUXls993hscegbduslSVJqp9M9zxoCyxI298WiMCLace+IgkZpEZLX6rR8ECSVBfjZ6RNltjJngdNKkb44x+rBgd77gmPP25wIEk5JtM9D74H1kvb3wOYBbyfdqwrrsCgDEnveVBSkr06JEm5w54HzSRG+P3v4ZprKo/tvTc88giUltbeTpLUImU6PBgJHB1COJWkB8J+wCMxxrK0a9YCvs3wc1WgHLYgSaqvCbOc86DJxQjnngvXXVd5bN994eGHDQ4kKUdletjCVcAc4EbgbyQBwtDykyGEnsCOwJgMP1cFyvBAklRfVYYtGB5kXoxw1llVg4MDDoARIwwOJCmHZbTnQYzxmxDChsDg1KEnYowT0i5ZA7gF+Gcmn6vCZXggSaqPGCPjZ1aGB2t0dthCRsUIZ5wBN99ceeygg+CBBxxfKEk5LtPDFogxTgb+Wsu5/wL/zfQzVbgMDyRJ9TFjwQzmLJoDQLuSdnRr6zKBGVM+x0F6cHDIIfCPfxgcSFIeyPSwBalZGR5IkuqjymSJndcghJDFavLMZZdVnRzx0EPhn/80OJCkPJHxngcAIYTewC4kSzLWNLgtxhgva4pnq7AYHkiS6iM9PHC+gwz6y19gyJDK/f33h/vug+Im+VVTkpQFGf8veghhGHBBtXsHIFbbNjxQoy1eXLlteCBJWpH0+Q4MDzLkttuSlRXK7b47PPigPQ4kKc9kdNhCCOGXwMXAaJJJEwNwD3AkcCdQBjwA7JzJ56pwpfc88HcUSdKKVB+2oEa65x44+eTK/R12gEcfdVUFScpDme55cBLwHfCLGOOS1DjCcTHGB4AHQgiPAv8G7s/wc1WgHLYgSaoPex5k0MMPw3HHVe5vuSU89RS0a5e9miRJTSbTEyb+DHg6xrgk7Vir8o0Y43PAc8B5GX6uCpThgSSpPqr0POhiz4MG+/e/4cgjoaws2d94Y3jmGejYMbt1SZKaTKbDgxJgatr+fKBztWs+BDbJ8HNVoAwPJEn14YSJGfDSS3DwwbAk9beiddeFF16Abi57KUn5LNPhwSSgd9r+BGDjatesAixBygDDA0lSXS1auohJsycBUBSKWKXjKlmuKAe9+Sbstx8sXJjs9+uXhAk9e2a3LklSk8t0ePAuydCFcv8Btg8hHBVCaB9C2Bs4OHWd1GiGB5Kkuvpu1nfE1OJPfTr2oaSVM+3Wy2efwd57w7x5yf6qqybBwSqGMJJUCDIdHjwFbBhC6Jfa/xMwExgOzAKeIFmB4aIMP1cFyvBAklRX42c4WWKDTZwIe+wBU1OjU7t3T4Yq9Ou3/HaSpLyR0dUWYozDSYKC8v1vQwhbAOcA/YFxwK0xxv9l8rkqXIsXV267VKMkaXlcprGBZsyAX/wCxqfCl3btkgkT11svq2VJkppXppdqXEaM8Rvg1KZ+jgqTPQ8kSXXlZIkNsGAB7L8//C/1d5/iYhgxArbaKrt1SZKaXaOHLYQQSkMIb4UQXgoh1Pq33xBC69Q1byzvOqk+DA8kSXU1fqbDFupl6VL45S/hlVcqj/3977DnntmrSZKUNZmY8+CXwEDgLzHGxbVdFGNcBFwDbJlqIzWa4YEkqa4ctlAPMcIpp8C//lV57Oqr4de/zl5NkqSsykR4cBDwdYzx6RVdGGN8FvgCOCQDz5UMDyRJdeawhXq49FK4447K/bPOgnPPzV49kqSsy0R4sBkwqh7XvwJsmoHnSoYHkqQ6iTFW7XnQxZ4HtbrjDhg6tHL/yCPh2mshhKyVJEnKvkyEBz2AH+px/Q9A9ww8VzI8kCTVyU/zfmL+kvkAdC7tTKfSTlmuqIV64gk4+eTK/d12g7vvhqJMr+4tSco1mfhJMB/oUI/rOwALMvBcqcpSjYYHkqTaOFliHYwdC0ccAWVlyf7AgfDII/6AlSQBmQkPvgW2qMf1mwMTVniVVAfpPQ9KXMNDklSLcTPGVWw7ZKEG334L++4L8+Yl+/36wdNPQ8eO2a1LktRiZCI8GAVsHULYfEUXhhAGAtsAIzPwXMlhC5KkOvli6hcV22t1XSuLlbRAs2fDPvvApEnJfpcu8O9/Q8+eWS1LktSyZCI8+CsQgYdDCOvXdlEIYT3gYWApcGsGnisZHkiS6uTzaZ9XbK/TfZ0sVtLCLFkChx0GH3yQ7BcXJ8szrl/rr3SSpAJV3NgbxBg/CyFcCgwF3g0hjAD+A3xHEiqsCuwCHAyUApfEGD9r7HOlGB22IEmqm/SeB2t3XzuLlbQgMcIZZ8Azz1Qeu/NO2Gmn7NUkSWqxGh0eAMQYLw0hLAGGAEcCR1S7JACLgQtjjFdl4pnS0qXJ7z0ArVolL0mSavL5VHseLOPGG+HWtM6gF14IxxyTtXIkSS1bRsIDgBjjlSGEfwDHAdsCvUlCg4nAq8DdMcbxy7mFVC8OWZAk1cWMBTP4cd6PALQpbsOqnVbNckUtwOOPw9lnV+4fdhhcemn26pEktXgZCw8AUuHAkEzeU6qN4YEkqS6qTJbYbS2KQiamfMphY8fCkUdWdt/bZhsYPhyKCvz7IklaLn9KKGctXly5bXggSarNF9PS5jvoVuDzHVRfknHNNeGxx6BNm6yWJUlq+QwPlLOcLFGSVBfOd5Aybx7sv/+ySzKutFJWy5Ik5QbDA+Ushy1IkurCngckQxSOOw7efTfZL1+Scb31sluXJClnGB4oZxkeSJLqIr3nQcEu0/inP8GDD1bu33KLSzJKkurF8EA5y/BAkrQiMcYqEyYW5LCFJ59MlmEsd/LJcMIJ2atHkpSTDA+UswwPJEkr8uO8H5m5cCYAHVp3oFf7XlmuqJl9/DH88peVKyvsuCPccENWS5Ik5SbDA+UswwNJ0opU73UQQshiNc1s+vRkgsTZs5P9NdaAhx92lmFJUoMYHihnGR5IklakynwHhTRZ4pIlcPjh8OWXyX67dvD4466sIElqMMMD5azFiyu3/SOKJKkm6SstFNR8B7//PTz/fOX+vffCJptkrx5JUs4zPFDOsueBJGlFCnKZxnvvheuuq9y/+GI4+ODs1SNJyguGB8pZhgeSpBVJH7ZQED0P3nyz6koKBxwAQ4dmqxpJUh4xPFDOMjyQJC1PWSzjy2lfVuyv3T3Pex788AMcdBAsXJjsb7hh0guhyF/3JEmN508T5SzDA0nS8kycPZF5i+cB0L1td7q17ZblippQ+QSJEycm+127JhMkduyY3bokSXnD8EA5y/BAkrQ86cs05n2vgwsvhFGjku0Q4P77oX//rJYkScovhgfKWYYHkqTlKZhlGh99FK6+unJ/6FDYY4+slSNJyk+GB8pZ6Us1Gh5IkqoriGUaP/8cjj66cn+vveCii7JXjyQpbxkeKGel9zwoKcleHZKklinvex7MnZsswTh7drLfty/8v//nBImSpCbhTxflLIctSJKWJ697HsSYLMn44YfJfmkpPPIIdMvjSSElSVlleKCcZXggSarNkrIlfDXtq4r9tbqtlcVqmsAtt8A//1m5f+utMGBA9uqRJOU9wwPlLMMDSVJtJsycwOKyZHKc3h1607E0j5YsHDMGzjqrcv+3v4XjjstePZKkgmB4oJxleCBJqk2V+Q7yaZnGH36AQw6BJUuS/YED4eabs1uTJKkgGB4oZxkeSJJq88XUtPkOuuXJfAdLl8IRR8DEicl+t24wYgS0aZPduiRJBcHwQDnL8ECSVJu87Hlw6aUwcmSyHQL84x/JCguSJDUDwwPlrMWLK7cNDyRJ6fJupYX//Acuu6xy/+KL4Re/yF49kqSCY3ignJXe86CkJHt1SJJanvTwYO1uOd7z4Icf4Je/TJZnBNhpJ7jkkuzWJEkqOAUVHoQQtg8hPBJCmBRCWJh6fz6EsFcN124TQng6hDAthDAvhPBBCOHMEEKr5dz/6BDCWyGEOSGEmSGEUSGEfZZzfdsQwrAQwmchhAUhhCkhhIdCCOsvp82qIYS7QggTU1/DuBDCDSGErvX/juQ2hy1IkmqyaOkixs0YB0Ag0L9b/+wW1BhlZfCrX8Hkycn+SislwxVa1frriCRJTaJgwoMQwkXAK8AOwLPAX4Anga7AoGrX7p927aPALUBr4HrggVrufy0wHOgN3AncB/wMeDKEcGoN15cCLwCXALOAG4EXgQOBt0MIW9XQpj8wFjgWeCtVz9fAGcDrIYTudftu5AfDA0lSTb6e/jVlsQyA1TuvTpviHJ5Q8Kqr4MUXk+0Q4L77oHfv7NYkSSpIxdkuoDmEEA4BLiP5cH5QjHF2tfMladudSD78LwUGxRjfTh2/GPgPMDiEcHiM8YG0NtsA5wBfAVvEGKenjl9D8mH/2hDCUzHGcWmPPRvYFhgBHBZj8ltOCOFB4DHgrhDCz8qPp9wK9AROjzFWrMsUQrgOOAu4AjixQd+kHGR4IEmqSd5MlvjKK1WHJ/zhD7D77tmrR5JU0PK+50EIoQj4MzAPOLJ6cAAQY0ybeo/BwErAA+XBQeqaBcBFqd2Tqt2i/AP7FeXBQarNOJJeC6UkvQXKawppbc5PDwhijI8Do4ENgB3T2qwJ7A6U3zPdEGAucFQIoX31ry9fGR5IkmqSF8s0/vhjsixjWepXhO22g2HDsluTJKmg5X14AGwD9AOeBqaHEPYOIfw+hHBGCOHnNVy/c+r92RrOvUISQmyTGnZQlzbPVLsGoD+wOvB5jPGbOrYp336+Wm8EUoHIa0A7YOsa7peXDA8kSTXJ+Z4HZWVw9NEwcWKy37073H8/FBdEh1FJUgtVCD+Ftki9/wC8QzIPQYUQwivA4Bjjj6lD66beP6eaGOOSEMI3wIbAmsAnqb/0rwLMiTFOquH55X/+SP/TR63PaGSb3VNtXqrlGgBCCGNrObXe8tq1NC7VKEmqSc4v0/iXv8Azz1Tu33MPrLpq9uqRJInC6HnQM/V+ItAW2BXoCGwEPEcyKeLDadd3Tr3PrOV+5ce7NPD65myT1+x5IEmqSZWeB7m2TOPrrydzG5Q791zYe+/s1SNJUkoh9DwoX8sokPQweD+1/1EI4UCSv+TvGEL4eYzx9TrcL6TeYz3rqM/1DXlGndvEGAfWeIOkR8KAejwzq9LDg5KS2q+TJBWOuYvm8v3s7wEoLiqmb5e+2S2oPqZPh8MPh6VLk/2tt4Yrr8xuTZIkpRRCz4PyCQy/TgsOAIgxzifpfQCwZeq9/C/4nalZp2rXrej6mnoM1PcZDW2T1+x5IEmq7stpX1Zs9+vSj5JWOZIuxwi/+x1MmJDsd+mSzHNgOi5JaiEKITz4LPU+o5bz5eFC22rXLzNIMoRQTDL54hLga4AY41zge6BDCKGmhZfL+0umz1VQ6zMy3CavGR5IkqrL2fkOhg+Hh9NGUd51F/Ttm61qJElaRiGEB6+QfNhfO4RQ00fMjVLv41Lv/0m9/6KGa3cgWdFgTIxxYdrx5bXZs9o1AF8BE4B1Qgj96thmZOp999TykxVCCB2BbYH5wBs13C8vGR5IkqrLyfkOvvgCTjutcv+EE+DAA7NXjyRJNcj78CDG+BPwIEl3/0vSz4UQdgP2IOnqX77M4gjgJ+DwEMLmade2AS5P7d5W7TG3p94vDCF0TWvTFzgFWAjcnVZTTGtzdXoYEELYH9ge+Bh4Oa3NV8DzQPk90w0D2gP3pnpCFATDA0lSdTnX82DRIjjySJib+vG97rpw3XXZrUmSpBoUwoSJAGcDW5F8uN8BeAtYAzgQWAocH2OcARBjnBVCOJ4kRBgVQngAmAbsR7Jc4giSMKJCjHFMCOG61HM+CCGMAFoDhwHdgNNijOOq1XQdsA8wGHgzhPASsDpwCDAPOC7GWFatzcnAGOCmEMIuwCepr2snkuEKFzb0G5SLDA8kSdVV6XnQPQd6HgwZAm+/nWyXlCTzHLRvn92aJEmqQd73PACIMU4h+ZB9PbAacDqwM/BvYPsY48PVrn8M2JFkyMPBwGnAYpJw4PBUz4HqzzgHOAaYDJwA/Br4CNg3xvjXGq5fSLJs5KUkyyueBewGPAZsEWN8s4Y2XwGbA8NTX885QH/gJuDnMcapdfyW5IXFiyu3DQ8kSQBfTK3sedDihy2MHAl//nPl/pVXwmabZa8eSZKWI9TwOVgFKoQwdsCAAQPGjh2b7VLqpHXrygBh4UIDBEkqdDMWzKDrn5PRg6WtSpl34TyKQgv9O8m0abDxxvB9sqwku+4Kzz0HRS20XklSXhg4cCDvvPPOOzHGgfVt608o5aQYq/Y8cCUrSVJ6r4O1uq3VcoODGJNJEcuDg+7d4Z57DA4kSS2aP6WUk6oHByFkrxZJUsuQM5Ml3nUXPPJI5f7f/w59+mSvHkmS6sDwQDnJyRIlSdXlxDKNn38Op59euX/iibD//tmrR5KkOjI8UE4yPJAkVdfiex6UL8s4b16yv/768Je/ZLcmSZLqyPBAOcnwQJJUXYtfpnHYMCiflLh1a/jnP6Fdu+zWJElSHRkeKCcZHkiS0sUYq0yY2OJ6HowZA3/6U+X+VVfBpptmrRxJkurL8EA5KX3CRMMDSdKP835k5sKZAHRo3YFe7XtluaI0c+bAr38NZWXJ/k47wZlnZrUkSZLqy/BAOSm954HLNEqSqvc6CC1pGZ7zz4evvkq2O3WC4cNdllGSlHP8yaWc5LAFSVK6FrvSwrPPwm23Ve7fdBOsvnr26pEkqYEMD5STDA8kSek+/enTiu0WEx5MmwbHHVe5f+CByfAFSZJykOGBcpLhgSQp3QdTPqjY/lmvn2WxkjQnnwyTJiXbPXvCHXdASxpOIUlSPRgeKCcZHkiS0n3wQ2V4sEmvTbJYScoDD8CDD1bu33knrLRS9uqRJKmRDA+UkwwPJEnlfpr3ExNnTwSgTXEb1uq2VnYL+v57OOmkyv3jjoP99stePZIkZYDhgXKSSzVKksr974f/VWxv1HMjWhW1yl4xMSZhwYwZyX7fvnD99dmrR5KkDDE8UE6y54Ekqdz7P7xfsb1xz42zWAlw++3w/PPJdgjJsoydOmW1JEmSMsHwQDkpPTwoKcleHZKk7Euf72DjXlkMD774As49t3L/7LNhxx2zV48kSRlkeKCcZM8DSVK5KpMlrpylyRKXLoVjjoF585L9DTeEyy/PTi2SJDUBwwPlJMMDSRLAkrIlfPTjRxX7P+uZpWUab7gBxoxJtouL4f/9P2jTJju1SJLUBAwPlJMMDyRJAF9O+5IFSxYAsErHVejernvzF/HZZ3DRRZX7F10Em23W/HVIktSEDA+UkwwPJEkA709OmywxG/MdLF0Kxx4LC5IAg003hT/+sfnrkCSpiRkeKCcZHkiSoAVMlnjDDfD668l2cXGyuoIz+UqS8pDhgXLS4sWV24YHklS4PpiSNllir2aeLLH6cIWLL4ZNsjRhoyRJTczwQDnJpRolSZDFngc1DVf4wx+a7/mSJDUzwwPlJIctSJJmLJjBhJkTAGjdqjXrdF+n+R5+/fUOV5AkFRTDA+UkwwNJUnqvgw1W2oCSVs304f3TTx2uIEkqOIYHykmGB5KkrAxZKB+usHBhsr/ZZg5XkCQVBMMD5STDA0lSenjQbJMlXncdvPFGsl1S4nAFSVLBMDxQTjI8kCQ1e8+DTz9NhiiUu/hi2DgLy0NKkpQFhgfKSYYHklTYymIZ/5vyv4r9Jg8Pqg9XGDAALrigaZ8pSVILYnignLR4ceW24YEkFZ6vpn3FvMXzAOjVvhc92/ds2gf+9a9VhyvcfbfDFSRJBcXwQDkpveeBv7tJUuFp1iEL33wDf/xj5f6FFzpcQZJUcAwPlJMctiBJha3ZJkuMEU44AeYlvRzYaCNXV5AkFSTDA+UkwwNJKmwfTGmmngd33w0vvphsFxXB3//uDx5JUkEyPFBOMjyQpMLWLMMWJk6Es8+u3D/rLNhyy6Z5liRJLZzhgXKS4YEkFa5ZC2fx9fSvASguKma9Hutl/iExwimnwMyZyX7//nDppZl/jiRJOcLwQDnJ8ECSCteHUz6s2F6vx3qUFpdm/iGPPAKPPVa5f+ed0K5d5p8jSVKOMDxQTnKpRkkqXE0+WeK0aUmvg3InnAA77ZT550iSlEMMD5ST7HkgSYWryec7OPtsmDIl2e7TB66+OvPPkCQpxxgeKCelhwclJdmrQ5LU/Jo0PHj2Wbjnnsr922+Hzp0z+wxJknKQ4YFykj0PJKkwlcWypgsPZs+G3/2ucv+II2DffTN3f0mScpjhgXKS4YEkFabxM8Yze9FsALq37U7vDr0zd/M//hEmTEi2u3eHG2/M3L0lScpxhgfKSYYHklSYqkyWuPImhBAyc+MxY+CWWyr3b7oJVlopM/eWJCkPGB4oJxkeSFJhqjJkoWeGhiwsWpSsqBBjsr/33smQBUmSVMHwQDnJ8ECSCtMHU5pgvoOrr4aPPkq227eHW2+FTPVokCQpTxgeKOcsXQplZcl2URG0apXdeiRJzef9ye9XbGckPPjsM7jsssr9K66A1Vdv/H0lScozhgfKOYsXV27b60CSCsfcRXP5ctqXABSFIjZYaYPG3bCsLBmuUN6dbYst4NRTG1mlJEn5yfBAOSd9yEJJSfbqkCQ1r49+/IhIMi/But3XpW1J28bd8K674JVXku1WreDOO+3OJklSLQwPlHOc70CSClOVyRIbO2Rh8mQ477zK/fPOg002adw9JUnKY4YHyjmGB5JUmDIaHpxxBsyYkWz37w+XXNK4+0mSlOcMD5RzDA8kqTC9/0OGJkt86il46KHK/TvugLaNHAIhSVKeMzxQzjE8kKTCE2PMTM+D2bPh5JMr948+GnbZpZHVSZKU/wwPlHNcbUGSCs93s75jxoIZAHRp04XVOq3WsBtdfDF8+22y3aMHXHttZgqUJCnPGR4o59jzQJIKT/VeByGE+t/krbfgppsq92+4IQkQJEnSChkeKOe4VKMkFZ4q4UHPBgxZWLwYjj8eYrLUI3vsAUcemaHqJEnKf4YHyjn2PJCkwvPeD+9VbDdovoPrr4cPUgFE27Zw223QkN4LkiQVKMMD5RzDA0kqPG9+92bF9ma9N6tf43HjYOjQyv1hw6Bfv4zUJUlSoTA8UM4xPJCkwjJp9iTGzxwPQNvitmzSa5O6N44RTjkF5s9P9jfZBM48M/NFSpKU5wwPlHMMDySpsLz+3esV25v32ZySVvWY8GbECHj66WQ7BLjjDifMkSSpAQwPlHMMDySpsLz+bWV48PNVf173hjNnwhlnVO6fdBJstVUGK5MkqXAYHijnLF5cuW14IEn5L73nwc9Xq0d4cOGFMGlSsr3yynDllRmuTJKkwmF4oJxjzwNJKhyLli7i7YlvV+zXuefBW2/BrbdW7t94I3TunOHqJEkqHIYHyjnp4YHDViUpv703+T0WLl0IQL8u/ejVodeKGy1ZAieckEyWCLDnnnDIIU1YpSRJ+c/wQDnHngeSVDje+O6Niu06D1m48UZ4//1ku21buOWWZLJESZLUYIYHyjmGB5JUOKrMd1CXIQvjx8Mll1TuDxkC/fo1QWWSJBUWwwPlHMMDSSoc9VppIUY49VSYNy/Z32gjOPvsJqxOkqTCYXignGN4IEmFYdLsSYyfOR6AtsVt2bjXxstv8Oij8NRTlft33OHkOJIkZYjhgXKO4YEkFYb0IQub99mcklbLCQJmzYLTT6/cP+EE2GabJqxOkqTCYnignLN4ceW24YEk5a96DVm45BL4/vtku2dP+NOfmrAySZIKj+GBco49DySpMFSZLHF5Ky288w7cfHPl/vXXQ9euTViZJEmFx/BAOSc9PHAoqyTlp0VLF/H2xLcr9mvtebB0KZx4IpSVJfu77gpHHNEMFUqSVFgMD5Rz7HkgSfnvvcnvsXDpQgD6delHrw69ar7wb3+D//432W7dGm65BUJopiolSSocBRkehBCOCiHE1Ou3tVyzTQjh6RDCtBDCvBDCByGEM0MIrZZz36NDCG+FEOaEEGaGEEaFEPZZzvVtQwjDQgifhRAWhBCmhBAeCiGsv5w2q4YQ7gohTAwhLAwhjAsh3BBCKJj+mYYHkpT/qsx3UNuQhcmT4Q9/qNz/wx9gnXWauDJJkgpTwYUHIYTVgJuBOcu5Zn/gFWAH4FHgFqA1cD3wQC1trgWGA72BO4H7gJ8BT4YQTq3h+lLgBeASYBZwI/AicCDwdghhqxra9AfGAscCb6Xq+Ro4A3g9hNB9RV9/PjA8kKT8V2W+g9qGLJxzDsycmWyvtRZccEEzVCZJUmEqqPAghBCAu4GpwO21XNOJ5MP/UmBQjPE3McbzgE2B14HBIYTDq7XZBjgH+ArYOMZ4VozxFGAgMA24NoTQt9qjzga2BUYAW8UYfx9jPBIYDLQD7gohVP/f51agJ3B6jPGAGOMFMcadSUKEdYEr6vs9yUWGB5KU/1YYHrz4Ivzzn5X7t94Kbdo0Q2WSJBWmggoPgNOBnUn+cj+3lmsGAysBD8QYK2ZqijEuAC5K7Z5Urc2JqfcrYozT09qMI+m1UJp6JlARYpS3OT/GWJbW5nFgNLABsGNamzWB3YHye6Ybkvp6jgohtK/l68obLtUoSflt4uyJTJg5AYC2xW3ZuNfGVS9YsABOPrly/4gjYLfdmrFCSZIKT8GEB6l5BP4E3BhjfGU5l+6cen+2hnOvAPOAbVLDDurS5plq1wD0B1YHPo8xflPHNuXbz6eHDQAxxtnAayQ9Frau4X55xZ4HkpTf3vjujYrtLVbZgpJW1ZbWufpq+OKLZLtTJ7juumasTpKkwlSc7QKaQwihGPh/wATgjyu4fN3U++fVT8QYl4QQvgE2BNYEPkn9pX8VYE6McVIN90v9dkP6DE61PqORbXZPtXmplmsACCGMreXUestr11IYHkhSfqsyWWL1IQtffAFXXlm5f+WVsPLKzVSZJEmFqyDCA5JJCTcDtosxzl/BtZ1T7zNrOV9+vEsDr2/ONnkpPTwoKan9OklSbkqf72DrVdM61MUIp5wCC5MlHNl8czjxRCRJUtPL+/AghLAlSW+Dv8QYX1/R9XW5Zeo91rNdfa5vyDPq3CbGOLDGGyQ9EgbU45lZYc8DScpfi5Yu4u2JFVMOVe158NBD8MILyXZREdx+O7SqdQVlSZKUQXk950HacIXPgYvr2Kz8L/idaznfqdp1K7q+ph4D9X1GQ9vkJcMDScpf701+j4VLk54F/br0o1eHXsmJmTPhzDMrLzzlFBhYYxYuSZKaQF6HB0AHkjkA1gcWhBBi+YtkhQKAO1PHbkjtf5Z6X4dqUmFEP2AJ8DVAjHEu8D3QIYTQu4Ya1k69p89VUOszMtwmLxkeSFL+qjLfwWppvQ4uuggmT062e/eGyy5r5sokSSps+T5sYSHw91rODSCZB+FVkg/m5b+t/Af4JfAL4P5qbXYgWdHglRjjwrTj/wGOSrW5u1qbPdOuKfcVyeSN64QQ+tWw4kJNbUam3ncPIRSlr7gQQugIbAvMB94gzxkeSFL+Sp/voGLIwtixcOutlRddfz10rq0jniRJagp53fMgxjg/xvjbml7AE6nL7kkdezC1PwL4CTg8hLB5+b1CCG2Ay1O7t1V71O2p9wtDCF3T2vQFTiEJMSpChRhjTGtzdQihKK3N/sD2wMfAy2ltvgKeB8rvmW4Y0B64N9UTIq8tXly5bXggSfllmfBg6dJkUsSyVGa+++5w6KFZqk6SpMKV7z0P6i3GOCuEcDxJiDAqhPAAMA3Yj2S5xBHAg9XajAkhXAecDXwQQhgBtAYOA7oBp8UYx1V71HXAPsBg4M0QwkvA6sAhwDzguPTeBSknA2OAm0IIuwCfAFsBO5EMV7iw8d+Bls+eB5KUnybOnsiEmRMAaFvclo17bQx33AFvpyZQLC2FW26BEJZzF0mS1BTyuudBQ8UYHwN2BF4BDgZOAxaThAOHp3oOVG9zDnAMMBk4Afg18BGwb4zxrzVcvxDYFbiUZHnFs4DdgMeALWKMb9bQ5itgc2A4SWhwDtAfuAn4eYxxakO/5lziUo2SlJ/S5zvYYpUtKPlxKvzxj5UX/OEPsNZaWahMkiQVbM+DGONQYOhyzr8G7FXPe94D3FOP6+eTTNw4ZEXXprX5Fji2PnXlG3seSFJ+WmbIwrnnJqssQBIa/P73WapMkiTZ80A5x/BAkvJTeniw73ft4R//qDx5663Qpk0WqpIkSWB4oBwTo8MWJCkfLVq6iLETxwLQeglsdeW9lScPOwx22y1LlUmSJDA8UI5ZsqRyu7gYivwXLEl54b3J77FwabIK8hXvdKX4iy+TEx07wnXXZbEySZIEhgfKMQ5ZkKT8VD5ZYr9pcPqLsypPXH459OmTpaokSVI5wwPllMWLK7cNDyQpf7z+3esQ4eZnoPWipcnBzTaDk0/ObmGSJAkwPFCOseeBJOWfGCOvffsaB34Ce3+ROhgC3H57MkZNkiRlneGBcoqTJUpS/vls6mdM//E7bnw27eDvfgdbbpm1miRJUlWGB8op9jyQpPzz3JfPMXQUrFY+1UHPnnDlldksSZIkVWNfQOUUwwNJyj+fjxrBjW+kHbj2WujaNWv1SJKkZdnzQDnF8ECS8svCRfP51W1jKI7J/vzttoJf/Sq7RUmSpGUYHiinGB5IUn755rpL+PmEMgAWtYK2f7s7mSxRkiS1KIYHyiku1ShJeeSnn1j9yr9W7P7n4AGw/vpZLEiSJNXG8EA5xZ4HkpRHzj+fdrMXAPBNF1hywe+zW48kSaqV4YFyiuGBJOWJV1+Fu++u2D1j7yJ22OAXWSxIkiQtj+GBckp6eFBSkr06JEmNsHgxnHRSxe4j68P0XbahU2mnLBYlSZKWx/BAOcWeB5KUB264AT78EIA5JXDmL2D3NXfPbk2SJGm5DA+UUwwPJCnHTZgAQ4dW7A7ZCb7rDLv3NzyQJKklMzxQTjE8kKQcd/rpMG8eAB/0hJu2gi5turB5n82zXJgkSVoewwPlFMMDScphTz4Jjz9esXvSPrCkFey65q60KmqVxcIkSdKKGB4opyxeXLlteCBJOWTuXDjttIrdZwatwpjVk+09+u+RpaIkSVJdGR4op9jzQJJy1OWXw/jxAMQePfjtNj9VnNptzd2yVZUkSaojwwPlFJdqlKQc9NFHcO21Fbsfn3s0E1svBGDd7uuyRpc1slWZJEmqI8MD5RR7HkhSjikrgxNPhCVLkv3ttuOeTSpPu8qCJEm5wfBAOcXwQJJyzPDh8OqryXZxMdx+O8+Ne6HitOGBJEm5wfBAOcXwQJJyyI8/wnnnVe6fdx6TVu/GBz98AEBJUQmD+g7KTm2SJKleDA+UUwwPJCmHnHceTJuWbPfrBxddxItfv1hxetvVt6VD6w5ZKk6SJNWH4YFyiks1SlKOGDUK7rmncv+WW6BdO57/+vmKQ7uv6ZAFSZJyheGBcoo9DyQpByxcmEySWO6QQ2DPPSmLZTz/VVp44HwHkiTlDMMD5RTDA0nKAddcA599lmx37Ag33ADABz98wJS5UwDo3rY7m/XeLEsFSpKk+jI8UE5JDw9KSrJXhySpFl9+CZdfXrl/xRXQpw9AlV4Hu/XfjaLgryGSJOUKf2orp9jzQJJasBjh5JOTYQsAm2+e7KdUGbLgfAeSJOUUwwPlFMMDSWrBHngAXngh2S4qgttvh1atAJi3eB6jJ4yuuNT5DiRJyi2GB8ophgeS1ELNmAFnnVW5f+qpMHBgxe4r419h0dLkP+IbrrQhq3RapZkLlCRJjWF4oJxieCBJLdQf/wg//JBs9+kDl11W5bSrLEiSlNsMD5RTFi+u3DY8kKQW4o03kiEK5W66CTp1qnLJc189V7FteCBJUu4xPFBOseeBJLUwixfD8ccnkyUC7LUXHHRQlUvGzRjHxz9+DEDrVq3ZYY0dmrtKSZLUSIYHyimGB5LUwvzlL/Dhh8l2u3Zwyy0QQpVLHvrooYrtXfrtQruSds1ZoSRJygDDA+WU9PCgpCR7dUiSgC+/hGHDKvcvuwz69l3msvTw4LAND2uGwiRJUqYZHiin2PNAklqIGOHEE2HBgmR/wAA4/fRlLvty2peMnTQWSIYs7L/e/s1ZpSRJyhDDA+UUwwNJaiHuuw9eeinZLiqCv/0NiouXuezhjx6u2N6j/x50adOlmQqUJEmZZHignGJ4IEktwE8/wVlnVe6fcQYMHFjjpQ9+9GDFtkMWJEnKXYYHyimGB5LUApx7LkydmmyvvjpcemmNl33202e8/8P7AJS2KmXfdfdtrgolSVKGGR4opyxeXLlteCBJWfDSS3DPPZX7t94KHTrUeGn6RIl7rb0XnUo7NXV1kiSpiRgeKKfY80CSsmj+/GSSxHKHHgp7713r5Q99XBkeHLrhoU1ZmSRJamKGB8ophgeSlEWXX54szwjQuTPceGOtl37848d8OOVDANoWt2WfdfZpjgolSVITMTxQzigrgyVLKvdbtcpeLZJUcD78EK6+unL/6qth5ZVrvTx9yMLe6+xNh9Y1D22QJEm5wfBAOaP6fAchZK8WSSooZWVwwgmVCe6228Jvf1vr5TFGV1mQJCnPGB4oZzhkQZKy5Lbb4PXXk+2SEvjb36Co9l8hPpzyIZ/+9CkA7Uvas9faezVHlZIkqQkZHihnGB5IUhZMmAAXXFC5f8EFsMEGy22S3utg33X3pV1Ju6aqTpIkNRPDA+UMl2mUpGYWY7K6wpw5yf5668GFF66gSawy38GhG7jKgiRJ+cDwQDnDngeS1Mz+8Q945plkOwT4+9+htHS5Td6b/B5fTPsCgA6tO7Dn2ns2dZWSJKkZGB4oZxgeSFIzmjIFzjijcv/UU2GbbVbYLL3Xwf7r7k+b4jZNUZ0kSWpmhgfKGenhQUlJ9uqQpIJwxhkwbVqyvfrqcOWVK2ziKguSJOUvwwPlDHseSFIzeeIJeOCByv2//Q06dFhhs7GTxvLNjG8A6Fzamd37795UFUqSpGZmeKCcYXggSc1g5kw46aTK/V//GvbYo05NH/ywstfBAesdQGnx8udHkCRJucPwQDnD8ECSmsH558PEicl2z55w3XV1ahZj5KGP01ZZ2NBVFiRJyieGB8oZhgeS1MRGjkyGKJT761+he/c6NX3z+zeZMHMCAF3bdGXXNXdtigolSVKWGB4oZyxeXLlteCBJGTZvHhx/fOX+AQfA4MF1bp6+ysKB6x1I61b+h1qSpHxieKCcYc8DSWpCQ4fC/2/vzuOsrurHj7/eM6wxAsoiaoJoimguYWpq7uY3lzQXwp9rluKauX2zb5pmfTPNXSk1yQy1UrHUcM0d11K+ihkIomigpIiyMzAz5/fH587MZZwZmJl7me317HH9fM7y+Zxz6Tzmzn3P+ZwzY0Z23qcP/OpXELFal1alqpWCB6O+6C4LkiR1NAYP1G4YPJCkInn5Zbjyytr0FVfA+uuv9uUT353I7IWzAejXsx97brRnoXsoSZJamcEDtRv5wYOuXVuvH5LUoZSXw7e/DVVVWXrPPeG7323SLW54+Yaa85FbjKRrqT+kJUnqaAweqN1w5oEkFcFPfwpvvJGd9+wJN9+82o8rAMxZNId7ptxTkz5l+1MaqS1JktorgwdqNwweSFKB/eMfcOmltenLLoNNNmnSLW5+5WYqqioA+Orgr7L1ulsXsoeSJKmNMHigdsPggSQV0LJlcNxxtY8r7L47nHZak25RUVXBTa/cVJM+9cunFrKHkiSpDTF4oHbD4IEkFdBFF8GUKdl5r15wyy1Q0rRfC+5/8/6ahRIH9hrIYVscVuheSpKkNsLggdqNFStqzw0eSFILvPBCtqNCtcsvh403bvJtfvWPX9WcnzjiRLqV+sNZkqSOyuCB2g1nHkhSASxduvLuCnvvDSed1OTbTPloCk+88wQAJVHCSds1/R6SJKn9MHigdsPggSQVwAUXwLRp2XlZGfz2t01+XAFW3p7x4GEHs2GfDQvVQ0mS1AYZPFC7kR886OoW4pLUdM8+C1dfXZu+6ioYMqTJt1m0fBG/f+33NelTt3ehREmSOjqDB2o3nHkgSS2wZAkcfzyklKX33RdOOKFZt7pj8h0sKF8AwLB+w9h76N6F6qUkSWqjDB6o3TB4IEkt8KMfwVtvZee9e8PYsRDR5NuklFZaKPGUL59CNOM+kiSpfTF4oHbD4IEkNdPTT8O119amr7kGNmzeGgXP/fs5Xv/wdQA+1/VzHLftcQXooCRJaus6fPAgIvpFxAkR8ZeIeCsilkbE/Ih4NiK+GxH1/htExM4R8WBEzIuIJRExOSLOjIjSRto6LiL+HhGLcm08FREHNlK/Z0RcHBFvRsSyiPgwIu6KiOGNXPP5iLglIt6PiPKImBkR10TE2k37l2l/3KpRkpphwYLscYVq+++f7bbQTPmzDo7e6mj69ujb/L5JkqR2o8MHD4CRwM3AjsBLwDXAPcAXgbHAXVFnvmVEHAw8A+wG/AX4FdANuBr4U32NRMQVwK3Aern2bge2Av4aEafXU7878DfgQmABcC3wGHAI8HJE7FjPNZsArwDHA3/P9edt4PvACxHRb7X+RdopZx5IUjOccQa880523rcv/OY3zXpcAWDOojnc8697atIulChJUufRpbU7sAZMAw4CHkgpVVVnRsSPyL6AHwYcShZQICJ6k335rwT2SCm9nMv/MfAEcHhEHJFS+lPevXYGzgFmANunlD7J5V9O9mX/ioiYkFKamdevs4FdgPHAqOq+RcSdwL3ALRGxVX6fgV8DA4EzUkrX57V/FXAW8HPg5Ob/U7VtBg8kqYnuvht+X7srAjfeCBts0OzbjZ00lhVV2TSwXTbchW0GbdPSHkqSpHaiw888SCk9kVL6a50v4aSU5gA35pJ75BUdDgwA/lQdOMjVXwZckEueUqeZ6i/sP68OHOSumUk2a6E72WwBAHIzHaqv+UF+31JK9wETgS2A3fOu2RjYF6i+Z76LgMXAMRHRq+6/QUfhVo2S1ASzZ8NJJ9Wmjz4aRo1q9u0qqiq46ZWbatLOOpAkqXPp8MGDVah+ir4iL2+v3PHheuo/AywBds49drA61zxUpw7AJsBgYFpK6Z3VvKb6/NF6AiELgeeAzwFfqed+HYIzDyRpNVVVZesafJKLZw8eDGPGtOiW9795P7MWzAJgYK+BHDb8sBZ2UpIktSed4bGFekVEF+DYXDL/S/+w3HFa3WtSShUR8Q6wJbAxMCX3l/4NgEUppQ/qaWp67rjZ6rTRwmv2zV3zeAN1AIiIVxoo2ryx61qbwQNJWk3XXguPPZadR8Btt0GfPi265a//8eua8xNHnEj3Lt0bqS1Jkjqazjzz4FKyRRMfTCk9kpdf/dvV/Aauq87v28z6a/KaDsXggSSthtdfhx/+sDZ93nmw224tuuXUuVN5/J0sLl0SJYzebnSL7idJktqfTjnzICLOIFvgcCpwTFMvzx1TE69rSv3mtLHa16SUtqv3BtmMhBFNaHONMnggSauwbBkcdVTtD8wRI+Dii1t8218+98ua84OGHcTgPoNbfE9JktS+dLqZBxFxGtm2iP8C9kwpzatTpfov+A3N7+xdp96q6tc3Y6CpbTT3mg5lxYrac4MHklSP88/PZh4A9OwJd9zR4h+Y0z6exrjXxtWkz/rKWS26nyRJap86VfAgIs4ExgD/JAsczKmn2pu542Z1C3LrJAwlW2DxbYCU0mJgNlAWEevVc79Nc8f8tQoabKPA13QozjyQpEY8/jhcdVVt+oorYPOWL2Vz8dMXU5kqAdhn433YbUjLHoGQJEntU6cJHkTEecDVwKtkgYMPG6j6RO749XrKdiPb0eD5lFL5al6zX506ADOA94DNImLoal7zZO64b0Ss9P9bRKwF7AIsBV6s534dgsEDSWrAvHlw3HG16f32g1Pq7ircdG98+AZ/fP2PNemf7fmzFt9TkiS1T50ieBARPyZbIPEVYO+U0txGqo8H5gJHRMSX8+7RA/jfXPKGOtfcmDueHxFr512zEXAaUA78rjo/pZTyrvllfjAgIg4GdiV7rOLpvGtmAI8C1ffMdzHQCxiXmwnRIeUHD7p2bb1+SFKbkhKcfDLMnp2l+/eHW27JdllooYueuoiUW0rngE0P4Cuf77C7AUuSpFXo8AsmRsRxwE+BSmAicEZ89heqmSmlWwFSSgsi4kSyIMJTEfEnYB5wENl2ieOBO/MvTik9HxFXAWcDkyNiPNANGAWsA3wvpTSzTptXAQcChwMvRcTjwGBgJLAE+E5KqarONacCzwPXRcTewBRgR2BPsscVzm/SP04748wDSarHb34Dd99dmx47FgYNavFtX53zKvdMuacm/dM9f9rie0qSpParwwcPyNYoACgFzmygztPArdWJlNK9EbE72Zfxw4AewFtkwYHrcjMHVpJSOiciJgOnA6OBKmAScHlKaUI99csjYh/gh8CRwFnAAuBe4KKU0r/quWZGbjbET8kekdgf+AC4Dri4nsUfOxSDB5JUx2uvwfe/X5sePRoOPrggt77wyQtrzg/Z/BBGrNdmN+ORJElrQNTzPVidVES8MmLEiBGvvPJKa3elXr17w8KF2fn8+VlakjqtRYtgu+1gWm6d3K22gpdeynZZaKGXZr3EV36bPaIQBK+d/BpbrbtVi+8rSZJa13bbbcekSZMmpZS2a+q1nWLNA3UMbtUoSTkpZQsiVgcOevWCu+4qSOAA4MKnamcdjPriKAMHkiTJ4IHaDxdMlKSc3/0Obr+9Nn3DDQXZlhFg4rsTeXTGowCURAk/2f0nBbmvJElq3wweqF2orISq3PKRpaXZS5I6pTfegNNPr00ffzwcc0xBbp1S4sdP/rgmfczWxzCs/7CC3FuSJLVvBg/ULrhYoiQBixfDyJGwdGmW3mILuP76gt3+iXee4Ol3s12Cu5R04cLdL1zFFZIkqbPoDLstqAPo2hUeeSQLIlTV3cBSkjqL00+HKVOy8549s3UOevUqyK1TSlzw5AU16eO3PZ6N1964IPeWJEntn8EDtQtdusC++7Z2LySpFY0bB7feWpseMwa23LJgt3/orYd4cdaLAHQr7cYFu12wiiskSVJn4mMLkiS1dVOnwqmn1qaPPjpb66BA6q51MHrEaAb3GVyw+0uSpPbP4IEkSW3ZkiXwrW9l6x0ADBuW7a4QUbAm7p16L5M+mARAjy49+NGuPyrYvSVJUsdg8ECSpLYqJTjhBHj99Szdo0e2zkFZWcGaWLJiCec8ek5N+rTtT2O9tdYr2P0lSVLHYPBAkqS26qqr4I9/rE1fdx1svXVBm7j4qYt559N3AFi7x9qct8t5Bb2/JEnqGAweSJLUFj32GPzgB7Xp0aPhxBML2sSrc17lyheurElfue+VDOg1oKBtSJKkjsHggSRJbc0778CoUbV70+60UzbroIAqqyo58a8nUpkqAdhjoz349rbfLmgbkiSp4zB4IElSW7JkCRxyCMybl6UHDYLx46F794I2M+bvY3j5/ZcB6F7anZsOvIko4CKMkiSpYzF4IElSW1G9QOJrr2Xprl3hnntg/fUL2sx789/j/CfOr0lfsNsFbNZvs4K2IUmSOhaDB5IktRV1F0gcMwZ23rmgTaSUOPWBU1m8Itv6cYsBW/CDXX6wiqskSVJnZ/BAkqS2oL4FEkePLngz4/81ngemP1CTvvkbN9OttFvB25EkSR2LwQNJklrbGlggEeCTpZ/wvYe+V5M+5cunsPOGhZ3ZIEmSOiaDB5IktabFi1deIHG99bJ1Dgq8QCLADx/7If9Z/J+smbL1+MXevyh4G5IkqWMyeCBJUmuprIQjj/zsAonrrVfwpia+O5HfTPpNTXrM/mPo06NPwduRJEkdk8EDSZJay9lnw/3316Z/9avskYUCK68oZ/SE2vUTDh52MIdsfkjB25EkSR2XwQNJklrDNdesvK7BuefCiScWpalLJl7C1LlTASjrVsb1+11PRBSlLUmS1DEZPJAkaU3785+zWQfVRo6Eyy4rSlOPv/04/zvxf2vSl+x1CRv22bAobUmSpI7L4IEkSWvSiy/CUUdBSll6551h3DgoKfxH8r/n/5sj7jmCqpTt4rDr4F05dftTC96OJEnq+AweSJK0psyYAQcdBMuWZekvfAHuuw969Ch4U+UV5Yy8eyRzl8wFYFDZIO48/E5KS0oL3pYkSer4DB5IkrQmfPwx7L8/fPRRlu7XDx56CPr3L0pzZz1yFi/NfgmA0ijlrsPvYr21Cr+LgyRJ6hwMHkiSVGzLlsE3vwnTpmXp7t2zXRa+8IWiNDfutXHc8PINNekr9r2CXYfsWpS2JElS52DwQJKkYqqqgm9/G559NktHwO23Z2sdFMGrc17lpAkn1aRHbTmK7+/4/aK0JUmSOg+DB5IkFUtKcM45cOedtXmXXw6HH16U5j5Z+gmH3XUYyyqyNRWG9x/O2IPGui2jJElqMYMHkiQVy49/DNdcU5s+7bSVt2gsoKpUxbH3Hsvbn7wNwFrd1uLPo/5MWbeyorQnSZI6F4MHkiQVw89/nr2qHXYYXHtt9thCEVwy8RImTJtQk/7dwb9j8/6bF6UtSZLU+Rg8kCSp0K6+Gi64oDZ9wAHwhz9AaXG2SXzkrUe48MkLa9L/vfN/c9gWhxWlLUmS1DkZPJAkqZBuvHHlRxP22QfGj4du3YrS3KtzXuWIe44gkQDYY6M9uGTvS4rSliRJ6rwMHkiSVCi//z2cckpt+qtfhXvvhR49itLcvz76F1+77Wt8uuxTANZfa33+dNif6FLSpSjtSZKkzsvggSRJhXDXXfCd79Smt98eHngAevUqSnPTP57O3uP2Zu6SuQD07dGXB458gHXL1i1Ke5IkqXMzeCBJUkvdfz8cdRRUVWXpbbaBhx+G3r2L0ty7n77L3uP2Zs6iOQCUdSvj4aMeZttB2xalPUmSJIMHkiS1xMMPw8iRUFGRpYcPh0cfhXXWKUpzsxfMZq9xe/HvBf8GoGeXnjxw5APs+Pkdi9KeJEkSGDyQJKn5xo+Hgw6C5cuz9CabwGOPwcCBRWnuw8Ufss9t+/D2J28D0K20G/cdcR+7DdmtKO1JkiRVM3ggSVJz3HwzjBoFK1Zk6cGD4fHHYf31i9LcvKXz+NptX2Pq3KkAdCnpwviR4/naJl8rSnuSJEn5DB5IktRUl10Go0fXrnEwbBhMnAhDhhSlufnL5vNft/8Xk/8zGYCSKOEPh/6Bbwz7RlHakyRJqsvggSRJqyslOO88+OEPa/NGjMgCB4MHF6XJeUvnsf8f9ufl918GIAhuPfhWRm45sijtSZIk1ceNoCVJWh2VlXDyyTB2bG3e7rtnOy0UaVeF6R9P58A/Hsi0j6fV5N144I0cs80xRWlPkiSpIQYPJElalfJyOProbIHEagcdBHfeCT16FKXJp2Y+xaF3Hsonyz6pybv269cyervRRWlPkiSpMQYPJElqzKJFcOih8Le/1eYdd1w2A6FLcT5Gb/m/WzhpwklUVGXbP/bs0pNxh4zj8C0OL0p7kiRJq2LwQJKkhsyeDYccAv/4R23emWfClVdCSeGXDaqsquR/Hv8fLn/+8pq8QWWDuP+I+9l+g+0L3p4kSdLqMnggSVJ9Xnghm3EwZ05t3s9+BuefDxEFb27R8kUc/eejue/N+2ryth20LfcfcT8b9tmw4O1JkiQ1hcEDSZLq+t3vssURly/P0qWlMGZMllcEsxbM4ht//Aavznm1Ju+gYQdxx6F3UNatrChtSpIkNYVbNUqSVG3FCjjjDPjOd2oDB/36ZesdFClw8Pjbj7PDzTusFDg4d6dz+fO3/mzgQJIktRnOPJAkCWDuXPjWt+DJJ2vztt4a7r0Xhg4teHNLVizhvL+dx5h/jKnJ61LShRsOuIETRpxQ8PYkSZJawuCBJEmTJ8PBB8PMmbV5hx0Gt94KZYX/6/+Ls17k2L8cy/R502vy+vXsx10j72KvoXsVvD1JkqSW8rEFSVLnNn487LTTyoGDn/0M7r674IGD5ZXLOf/x89nlll1WChwcNOwg3jj1DQMHkiSpzXLmgSSpc1q0CM46C8aOrc0rK4Pbb89mIRTY5P9M5ti/HMtr/3mtJm+tbmtx3X7Xcdw2xxFF2MFBkiSpUAweSJI6nxdegGOOgRkzavM22QTuuw+23LKgTVVUVXDl81dy4VMXsrxyeU3+Hhvtwa0H38qQvkMK2p4kSVIxGDyQJHUeK1ZkjyT8/OdQVVWbf/jhcNNNsM46BWsqpcSD0x/kv//230yZO6Umv0eXHly696V8b8fvURI+PShJktoHgweSpM7hzTfh6KPh5Zdr83r3hjFjsvwCPjbw6pxXOefRc3jinSdWyt9+/e0Zd8g4Nu+/ecHakiRJWhMMHkiSOraU4IYb4NxzYenS2vzdd4ff/x6GFO6xgVkLZnHBExcw7rVxJFJNflm3Ms7f9XzO3flcupT40StJktoff4ORJHVc77wDp54KDz9cm9e1a/bYwtlnQ2lpQZpZWL6Qy567jKteuIqlFbUBitIoZfR2o7lo94tYt2zdgrQlSZLUGgweSJI6niVL4NJL4Ze/hPLy2vwtt4Q77oBttilIM/OXzWfspLH88vlf8uHiD1cq+8Zm3+CyfS5j+IDhBWlLkiSpNRk8kCR1HCnB+PFwzjnw73+vXHb22dmMgx49WtzMe/Pf49oXr+XmSTezcPnClcq+NOhLXLHvFew1dK8WtyNJktRWGDyQJHUM//wnnHEGPPnkyvnbbQfXXw877dTiJl55/xWufOFK7nrjLipT5Upln+/9eS7Z6xKO2vood1GQJEkdjsEDSVL79skncNFF8OtfQ2XeF/oBA+AXv4Djj4eS5n+Zr0pVPDT9Ia544QqemvnUZ8o37785Z3/lbI7Z5hh6dGn5rAZJkqS2yOCBJKl9WrgQbrwxW9dg7tza/NJSOP10+MlPoG/fZt9++sfTuW3ybdw2+TZmfjrzM+V7brQn5+x0Dvttup8zDSRJUodn8ECS1L58+mn2GMI118C8eSuX7bUXXHddtjBiM8xbOo87/3kn4yaP48VZL36mvDRKGfXFUZyz0zmMWG9Es9qQJElqjwweSJLah48+gquvhjFjslkH+YYMgSuvhEMPhYgm3ba8opwHpz/IbZNvY8K0CayoWvGZOn179OW7X/ouZ+x4BoP7DG7Ju5AkSWqXDB5Iktq299+HK66Am27KtmDMN3Qo/M//wLHHQvfuq33LDxd/yEPTH2LC9Ak88tYjn9kxAaBLSRcO2PQAjt3mWA7Y9AC6d1n9+0uSJHU0Bg8kSW1PSvD00/Db38Ldd0N5+crlm28O558PRxwBXVb9UZZSYvJ/JjNh2gQmTJ/AS7NeIpHqrbvDBjtw7NbHMuqLo+j/uf6FeDeSJEntnsEDSVLb8cEHcOutcMst8NZbny3fZhu44AI45JBsYcRGvDf/PZ559xmeefcZHnrrIWYtmNVg3aF9h3LkVkdyzNbHMKz/sBa+CUmSpI7H4IEkqXVVVMCDD8LYsdkxf7vFajvumAUNDjig3jUNUkq8+fGbTHx3Is+89wwT353Iu/PfbbDJkihhlw134cDNDuTAzQ5keP/hRBPXSpAkSepMDB5Ikta8FSvgmWfgr3+Fu+7KZhzU1acPHHkknHACjFh5Z4MPF3/IpA8mMemDSbz8/ss8+96zfLTko0ab7NujL/t9YT8O3OxAvv6Fr7NOz3UK+Y4kSZI6NIMHkqQ149NP4aGH4P77s+P8+fXX2313+O534bDDSD17MmvBLP7vzftrggWTPpjE7IWzV9lcr6692GnDndht8G7svtHu7LzhznQp8WNPkiSpOfwtSpJUHCnBm2/CI49kAYNnnskeUaiv6qBBzPvWN5i037a83Gs+Uz9+jCl3XM/UuVPr3QmhPuv0XIevDv4quw3ejV2H7MqXBn2JrqVdC/mOJEmSOi2DB5Kkwli2DF5+GZ57Lns9/zx8/HGD1ecNWIvnvtSP+zat5A/9P2Bp3AwvrV5TPbv0ZNtB2zJivRGMWG8EO2ywA1sM2IKSKCnQm5EkSVI+gwftUER8Hvgp8HWgH/ABcC9wcUrpk1bsmqTOoqICZsyAf/6T9OKLVD77DCWT/o+S5Ssavezv68Nfh8H9w2DyugshVj2roG+Pvmyz7jY1gYIR641gWL9hlJY0vtuCJEmSCsfgQTsTEZsAzwMDgfuAqcAOwPeBr0fELimlhv/UJ0mrkFJi0fJFfLLsE+YtnsvCaf+k6vXJlE6ZStm0mazz9hwGzfqEbiuqAAga/jCZ2xOeGwwPbAoTNoMPejfc7oa9N2Tz/pszvP9whg8YXnM+sNdAd0KQJElqZQYP2p9fkwUOzkgpXV+dGRFXAWcBPwdObqW+SWpllVWVLF6xmMXLF7No+SIWLV/EgvIFNa/55fNrz5fNZ9HSTyl5fw49Zv+Hsjkfs86HCxnw0RI2/DQx5FPYbD58rv5lCuo1bR14dnAWMHhuQ3izP1l0AejRpQeb992IoX2HZq+1h7JRLj2s/zDKupUV459EkiRJBWDwoB2JiI2BfYGZwK/qFF8EjAaOiYhzUkqL13D3JJF9eV9RtYIVlStYUbWC5ZXLWV65nPKK8przmrzKcsoryllWsYzyytyxTnrpiqUsrVjK0hVLWVKxJDuuWFKTVx0oWFK+iLR4Ed0XldO7nJpX32UwYAkMWFx7HL4EBi7OzvsthdLUvPc6ey3450B4bV14ZaNuzNxyfXquP4QNem/A+mXrc3LvDdhgrQ3YoPcGDO07lHXL1nVNAkmSpHbK4EH7slfu+GhKqSq/IKW0MCKeIwsufAV4fE13rpiqKit4/a+/bfF9Ulr1t6REajTd1DY+c79VlOcn88vyr6uvT/XVrc77TJt55TV1agtXui6RagoTqSav5g7V6bRy+4lEFVVkl+fq5dWt+V9+XvV1qSrrR1665pjLryI7r0yVWTqvTlVVFSlVUZmyYxVVVFVWZsdceVWurLKqkqqqrKzmPGXllZUVWTrvVZnqpKsqqKqoIFVVUllZQaTsD+0lKXuVVmVfzhs6dqmCbpXQtRK6VmXHbpXQvQrKKqF7JfRcAT0rao898s4/tyILEqxVDsX6Wv7pWl2ZPbgvc4cOYuFmQ6gYPoySLbdi7fU3ZqNeA9mxbF36dO/jowWSJEkdmMGD9mVY7jitgfLpZMGDzWgkeBARrzRQtHnzu1ZclSuWs80hPo0hFcOy/n1Z8fn1qBo8mNKNhtJ9483ouvEXYMgQGDKEvn360Le1OylJkqRWZfCgfemTO85voLw6v2/xuyKprUm9ekHv3kTv3pD/GjCg9jVw4Mrn/fvTo2tXerR25yVJktSmGTzoWKrnDDc6zz6ltF29F2czEkYUulOFECUlvLFxoRZTa/nU6qbfIZqSbFJrDZfW30g0VP6ZnKjvsFLZZ+7VQBs1/41GyoBs1nteSxE15flldfMaShNRM5U+CCjJymr+G0CU5JXl/ldSXSuIkhIigogSKCmhpKSkJp8ooaS0lCjJXiUlpVl+SUn2ZkpKoLS08VfXrrWvbt0+e969O/TsCT16ZMfqV3X6c5/LAgRlZUQXf6RLkiSpOPxNs32pnlnQp4Hy3nXqdRhduvVgyxmr3g9ekiRJklR4LnvdvryZO27WQPmmuWNDayJIkiRJktRksTqrz6ttiIhNgLfItmrcJH/HhYhYC/iALCA0oDlbNUbExz179lxn+PDhBeqxJEmSJKmtmDJlCkuXLp2XUurX1Gt9bKEdSSnNiIhHyXZUOA24Pq/4YqAXcFNzAgc5C5YuXcqkSZNmtqynRVO9G8TUVu2F5FhU2+A4VFvgOFRb4VhUW9AexuFGwILmXOjMg3YmN/vgeWAgcB8wBdgR2JPscYWdU0oft14Pi6d6i8mGFnyU1hTHotoCx6HaAseh2grHotqCjj4OXfOgnUkpzQC+DNxKFjQ4B9gEuA7YqaMGDiRJkiRJrcfHFtqhlNK/geNbux+SJEmSpM7BmQeSJEmSJKlRBg8kSZIkSVKjDB5IkiRJkqRGuduCJEmSJElqlDMPJEmSJElSowweSJIkSZKkRhk8kCRJkiRJjTJ4IEmSJEmSGmXwQJIkSZIkNcrggSRJkiRJapTBA0mSJEmS1CiDB2rzIuLzEXFLRLwfEeURMTMiromItVu7b+pYIqJfRJwQEX+JiLciYmlEzI+IZyPiuxFR78/MiNg5Ih6MiHkRsSQiJkfEmRFRuqbfgzqmiDgmIlLudUIDdRyHKoqI2DUi7omID3Kfwx9ExKMRsX89dR2HKriIOCA35mblPpvfjoi7I2KnBuo7DtUsEXF4RFwfERMjYkHuc/f2VVzT5PEWEcdFxN8jYlHud82nIuLAwr+jwoqUUmv3QWpQRGwCPA8MBO4DpgI7AHsCbwK7pJQ+br0eqiOJiJOBG4APgCeB94B1gUOBPsA9wMiU94MzIg7O5S8D7gTmAd8AhgHjU0oj1+R7UMcTERsCrwOlQBlwYkppbJ06jkMVRURcAPwMmAtMIPv52B/4EvBkSukHeXUdhyq4iLgM+AHwMXAv2Vj8AnAQ0AU4NqV0e159x6GaLSJeBbYBFgGzgM2BO1JKRzdQv8njLSKuAM7J3X880A04AlgH+F5KaUxh31XhGDxQmxYRjwD7AmeklK7Py78KOAu4KaV0cmv1Tx1LROwF9AIeSClV5eUPAv4ObAgcnlK6J5ffG3iLLLCwS0rp5Vx+D+AJYCfg/6WU/rRG34g6jIgI4G/AUODPwLnUCR44DlUsETESuAt4DDg0pbSwTnnXlNKK3LnjUAWX+/ydDXwEbJ1S+jCvbE+ysfVOSmnjXJ7jUC2SG1ezyMbR7mR/TKo3eNCc8RYROwPPATOA7VNKn+TyNwJeIfs9dPOU0swivcUW8bEFtVkRsTFZ4GAm8Ks6xRcBi4FjIqLXGu6aOqiU0hMppb/mBw5y+XOAG3PJPfKKDgcGAH+q/sDI1V8GXJBLnlK8HqsTOAPYCzie7GdefRyHKrjcY1qXAUuAI+sGDgCqAwc5jkMVwxCy7ysv5QcOAFJKTwILycZdNcehWiSl9GRKaXr+LNNGNGe8Vf/R8+fVgYPcNTPJvu90J/vMb5MMHqgt2yt3fLSeL3MLyaJ2nwO+sqY7pk6p+pfkiry86jH6cD31nyH7pXvniOhezI6pY4qI4cClwLUppWcaqeo4VDHsTDbj5UHgk9wz5+dFxPcbeM7ccahimA4sB3aIiP75BRGxG7AW2cyYao5DrUnNGW+NXfNQnTptjsEDtWXDcsdpDZRPzx03WwN9UScWEV2AY3PJ/B/2DY7RlFIF8A7Z85gbF7WD6nByY+42snU3frSK6o5DFcP2ueN/gElk6x1cClwDPB8RT0dE/l98HYcquJTSPOA8svWH/hURv4mIX0TEXcCjZI91nZR3ieNQa1KTxltutvQGwKKU0gf13K/Nf7fp0todkBrRJ3ec30B5dX7f4ndFndylwBeBB1NKj+TlO0ZVLBeSLUj31ZTS0lXUdRyqGAbmjieT/QK8D/AS2TTyK4H/Au6m9lEux6GKIqV0TUTMBG4BTswregu4tc7jDI5DrUlNHW/tfnw680DtWeSOrvqpoomIM8hWxJ0KHNPUy3NHx6hWW0TsQDbb4MqU0guFuGXu6DhUU1RvMRZkC8U+nlJalFJ6AziEbEGx3RvaKq8ejkM1S0T8gGxF+luBTcgWlNsOeBu4IyJ+2ZTb5Y6OQ60JzR1vbXZ8GjxQW1YdfevTQHnvOvWkgoqI04BrgX8Be+amT+ZzjKqg8h5XmAb8eDUvcxyqGKoX8no7pfRafkFuNkz1LKwdckfHoQouIvYgW7jz/pTS2Smlt1NKS1JKk8iCWLOBc3KLbIPjUGtWU8fbquqvamZCqzN4oLbszdyxoed+Ns0dG1oTQWq2iDgTGAP8kyxwMKeeag2O0dyXwKFkCyy+XaRuquMpIxtPw4FlEZGqX2S7zADcnMu7Jpd2HKoYqsfVpw2UVwcXetap7zhUIR2YOz5ZtyCltIRsG+USsse8wHGoNatJ4y2ltJgs4FUWEevVc782/93G4IHasuoPin1zW0bViIi1gF2ApcCLa7pj6tgi4jzgauBVssDBhw1UfSJ3/Ho9ZbuR7QbyfEqpvOCdVEdVDvy2gdf/5eo8m0tXP9LgOFQxPEP2S++mEdGtnvIv5o4zc0fHoYqhepX6AQ2UV+cvzx0dh1qTmjPeGrtmvzp12hyDB2qzUkozyFbS3Qg4rU7xxWTPvI3LRfGkgoiIH5MtkPgKsHdKaW4j1ccDc4EjIuLLeffoAfxvLnlDsfqqjieltDSldEJ9L+D+XLXf5/LuzKUdhyq43M++O8mm0V6YXxYRXyNbMHE+tTvQOA5VDBNzx9ERsUF+QUTsR/aHpGXA87lsx6HWpOaMtxtzx/MjYu28azYi+75TDvyuWB1uqUipza7HIBERm5B9IAwE7gOmADsCe5JN6dk5pfRx6/VQHUlEHEe2IFMlcD31P3M2M6V0a9413yT78FgG/AmYBxxEtn3PeOBbyR+0KoCI+AnZowsnppTG1in7Jo5DFVhEDASeA75A9iXu72S7LRxCtqDXkSmlu/PqfxPHoQooN/P0EbLdPhYCfwHmkD3adSDZgnRnppSuzbvmmzgO1Uy58fPNXHIQWaD0bWoDWXNTSufWqd+k8RYRVwJnky08Ox7oBowC+gHfSymNKfw7KwyDB2rzImJD4Kdk03v6AR8A9wIX17OAndRseV/OGvN0SmmPOtftApwP7AT0INs+6hbgupRSZeF7qs6oseBBrtxxqIKLiHWAC8gCBhuQfYF7FvhFSukzjw06DlVoEdGV7C+yRwBbkE0Fn0cWzLoupfRoPdc4DtUsq/G74LsppY3qXNPk8Zb7g9XpZGO6CpgEXJ5SmtDCt1BUBg8kSZIkSVKjXPNAkiRJkiQ1yuCBJEmSJElqlMEDSZIkSZLUKIMHkiRJkiSpUQYPJEmSJElSowweSJIkSZKkRhk8kCRJkiRJjTJ4IEmSJEmSGmXwQJIkSZIkNcrggSRJkiRJapTBA0mSJEmS1CiDB5IkqUOJiG9HRIqIb7d2XyRJ6ii6tHYHJEmSGhIRqYmXHF+UjkiS1MkZPJAkSW3ZxfXknQn0Aa4FPq1T9irwDvAi8EER+yVJUqcSKTU1oC9JktR6ImImMAQYmlKa2bq9kSSpc3DNA0mS1KE0tOZBRMzMvcoi4uqI+HdELI2IVyPim7k6XSLiRxExPSKWRcSMiDi9kbb+KyIejIi5EVGeq395RPQt6puUJGkN87EFSZLUmXQF/gasA9wHdAP+H3BPROwLnArsCDwElAMjgesj4qOU0p35N4qIC8keq5gHTAA+BLYGzgX2j4idUkoL1si7kiSpyAweSJKkzmR9YBKwR0qpHCAibgOeAe4GZgBfTCl9miu7CpgK/BCoCR5ExJ5kgYMXgP2r6+fKvg38Lld+VrHfkCRJa4KPLUiSpM7mzOrAAUBKaSLZIotrA+flBwJSSm8DzwFbRURp3j3OyB1PzK+fu+ZWsoUbjypC3yVJahXOPJAkSZ3JpymlGfXkvw8MBV6pp2w2UAoMyp0D7ASsAEZGxMh6rukGDIiIfimlj1vebUmSWpfBA0mS1JnMbyC/AiClVF95Re7YNS+vH9nvURetor0ywOCBJKndM3ggSZLUdPOBkpTSOq3dEUmS1gTXPJAkSWq6F4G1I2LL1u6IJElrgsEDSZKkprs6d7w5ItavWxgRvSLiK2u4T5IkFY2PLUiSJDVRSunxiPgh8AtgekQ8SLZjQxkwBNgdeBb4euv1UpKkwjF4IEmS1Awppcsi4jmybRu/ChxMthbCbOA3wB9asXuSJBVUpJRauw+SJEmSJKkNc80DSZIkSZLUKIMHkiRJkiSpUQYPJEmSJElSowweSJIkSZKkRhk8kCRJkiRJjTJ4IEmSJEmSGmXwQJIkSZIkNcrggSRJkiRJapTBA0mSJEmS1CiDB5IkSZIkqVEGDyRJkiRJUqMMHkiSJEmSpEYZPJAkSZIkSY0yeCBJkiRJkhpl8ECSJEmSJDXK4IEkSZIkSWqUwQNJkiRJktSo/w/UdIS7YMxqhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 370,
       "width": 519
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize= (8,6))\n",
    "tt= np.linspace(0,100,100)\n",
    "plt.plot(tt, logistic_model(tt,4.8,50,100000), c='b', label='Logisitc')\n",
    "plt.plot(tt, expo_model(tt, 0.23), c='g', label='Exponential')\n",
    "plt.plot(tt, gompertz_model(tt, 12, 49, 100000), c='r', label='Gompertz')\n",
    "plt.legend()\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Cases')\n",
    "plt.ylim(-1000,150000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### COVID-19 Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = '/Users/genie/data/covid-19/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/genie/data/covid-19/.DS_Store\n",
      "/Users/genie/data/covid-19/jhu_time_series_covid19_deaths_global_narrow.csv\n",
      "/Users/genie/data/covid-19/owid-covid-data.csv\n",
      "/Users/genie/data/covid-19/eiu_democracy_indices.csv\n",
      "/Users/genie/data/covid-19/jhu_time_series_covid19_confirmed_global_narrow.csv\n",
      "/Users/genie/data/covid-19/jhu_time_series_covid19_recovered_global_narrow.csv\n",
      "/Users/genie/data/covid-19/us_daily.csv\n"
     ]
    }
   ],
   "source": [
    "for dirname, _, filenames in os.walk(data_dir):\n",
    "    for filename in filenames:\n",
    "        print(os.path.join(dirname, filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 393 ms, sys: 43.5 ms, total: 437 ms\n",
      "Wall time: 465 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "data = pd.read_csv(os.path.join(data_dir,'owid-covid-data.csv'))\n",
    "data = data[data.iso_code.notnull()][['iso_code', 'location', \n",
    "       'date', 'total_cases', 'new_cases','total_deaths', 'new_deaths', 'total_tests','new_tests', 'population']]\n",
    "data = data.rename(columns={'location':'country_name'})\n",
    "\n",
    "# Merge recoveries data from JHU dataset\n",
    "recovered_df = pd.read_csv(os.path.join(data_dir,'jhu_time_series_covid19_recovered_global_narrow.csv'))\n",
    "recovered_df = recovered_df.rename(columns={'iso_country_code':'iso_code','recovered':'total_recoveries'})\n",
    "recovered_df = recovered_df.groupby(['iso_code','date'])['total_recoveries'].sum().reset_index()\n",
    "recovered_df['date'] = pd.to_datetime(recovered_df['date'])\n",
    "recovered_df['date'] = recovered_df['date'].dt.strftime('%Y-%m-%d')\n",
    "recovered_df['new_recoveries'] = recovered_df['total_recoveries'] - recovered_df.sort_values(by=['date'], ascending=True)\\\n",
    ".groupby(['iso_code'])['total_recoveries'].shift(1)\n",
    "data = pd.merge(data, recovered_df, how='left', on=['iso_code','date'])\n",
    "del recovered_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EXAMPLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# iso_code = 'IND'\n",
    "# cases_starting_point = 100\n",
    "\n",
    "# df = data[(data.iso_code == iso_code) & (data.total_cases >= cases_starting_point)][['date','total_cases','new_cases']]\n",
    "# df['time'] = df.reset_index().index\n",
    "\n",
    "# # cutoff_date = '2020-05-31'\n",
    "# # df = df[df.date <= cutoff_date]\n",
    "\n",
    "# ##remove last data point to make prediction\n",
    "# x = np.array(list(df.time.values))[0:-1]  \n",
    "# y = np.array(list(df.total_cases.values))[0:-1]\n",
    "\n",
    "# y0 = y[0]\n",
    "# yf = y[-1]\n",
    "\n",
    "# country_name = data[data.iso_code == iso_code].country_name.values[0]\n",
    "\n",
    "# print(\"~~ Prediction for {}-{} as of {} ~~\".format(country_name, iso_code, df.date.max()))\n",
    "# print(\"Initial no. of cases on {} -> {}\".format(df.date.min(), y0))\n",
    "# print(\"No. of cases on {} -> {}\".format(df.date.values[-2], yf))\n",
    "# print('>>> Exponential Model')\n",
    "# fit_exp = curve_fit((lambda tt,pp: expo_model(tt,pp,N0=y0)),x,y,p0=[0.10])  ##Take N0 from initial number of reported cases\n",
    "# pi, sigma_pi = fit_exp[0][0], np.sqrt(fit_exp[1]).flatten()[0]\n",
    "\n",
    "# print([ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('p','sigma_p'),(pi,sigma_pi))])\n",
    "\n",
    "# y_pred = expo_model(x,pi)\n",
    "# MSLE=sklm.mean_squared_log_error(y,y_pred)\n",
    "# print(\"Mean squared log error (MSLE): \", '{:.3f}'.format(MSLE))\n",
    "# print(\"Exp of RMSLE: \", '{:.3f}'.format(np.exp(np.sqrt(MSLE))))\n",
    "# print(\"R2 score: \", '{:.3f}'.format(sklm.r2_score(y,y_pred)))\n",
    "\n",
    "# ## Logistic\n",
    "# # print('\\n')\n",
    "# print('>>> Logistic Model')\n",
    "\n",
    "# # fit_i = curve_fit(logistic_model,x,y,p0=[3,50,5000], maxfev=10000)#, bounds=([0,0,0],[10,100,150000]))\n",
    "\n",
    "# fit_i = curve_fit(logistic_model,x,y,p0=np.random.exponential(size=3))#, bounds=([0,0,0],[10,100,150000]))\n",
    "\n",
    "# # fit_i = curve_fit(logistic_model,x,y,p0=[3,50,5000], bounds=([0,0,0],[10,100,150000]))\n",
    "\n",
    "# # fit_i = curve_fit(logistic_model,x,y, p0=np.random.exponential(size=3) , bounds=([0,0,0],[10,10,10000000000]))\n",
    "# ai,bi,ci = fit_i[0]\n",
    "# sigma_ai, sigma_bi, sigma_ci = np.sqrt(np.diag(fit_i[1]))\n",
    "\n",
    "# print([ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('a','b','c'),(ai,bi,ci))])\n",
    "# print([ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('sigma_a','sigma_b','sigma_c'),(sigma_ai,sigma_bi,sigma_ci))])\n",
    "\n",
    "# y_pred = logistic_model(x,ai,bi,ci)\n",
    "# MSLE=sklm.mean_squared_log_error(y,y_pred)\n",
    "# print(\"Mean squared log error (MSLE): \", '{:.3f}'.format(MSLE))\n",
    "# print(\"Exp of RMSLE: \", '{:.3f}'.format(np.exp(np.sqrt(MSLE))))\n",
    "# print(\"R2 score: \", '{:.3f}'.format(sklm.r2_score(y,y_pred)))\n",
    "\n",
    "# perc_flat = 0.98\n",
    "# sol = int(fsolve(lambda x : logistic_model(x,ai,bi,ci) - perc_flat*int(ci), bi))\n",
    "# datesol = datetime.strftime(datetime.strptime(df.date.min(),'%Y-%m-%d') + timedelta(days=sol), ' %d, %b %Y' )\n",
    "# print('Day of flattening of the infection curve -> {}'.format(datesol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** Prediction for China-CHN as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-01-19 -> 216.0\n",
      "No. of cases on 2020-07-28 -> 86783.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.035', 'sigma_p=0.001']\n",
      "MSLE: 45.617, Exp of RMSLE: 9.670, R2 score: -9.730\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.113', 'c=84166.683'] ['sigma_a=1.115', 'sigma_b=0.005', 'sigma_c=390.112']\n",
      "MSLE: 0.418, Exp of RMSLE: 1.909, R2 score: 0.963\n",
      "Day of flattening of the infection curve ->  13, Mar 2020\n",
      "\n",
      "\n",
      "**** Prediction for Thailand-THA as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-16 -> 114.0\n",
      "No. of cases on 2020-07-28 -> 3297.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.029', 'sigma_p=0.001']\n",
      "MSLE: 22.803, Exp of RMSLE: 9.670, R2 score: -11.185\n",
      "~~ Logistic Model ~~\n",
      "['a=9.005', 'b=0.139', 'c=3111.112'] ['sigma_a=0.666', 'sigma_b=0.004', 'sigma_c=10.831']\n",
      "MSLE: 0.030, Exp of RMSLE: 1.190, R2 score: 0.984\n",
      "Day of flattening of the infection curve ->  28, Apr 2020\n",
      "\n",
      "\n",
      "**** Prediction for South Korea-KOR as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-02-21 -> 155.0\n",
      "No. of cases on 2020-07-28 -> 14203.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.032', 'sigma_p=0.001']\n",
      "MSLE: 31.421, Exp of RMSLE: 9.670, R2 score: -11.591\n",
      "~~ Logistic Model ~~\n",
      "['a=3.858', 'b=0.074', 'c=12089.302'] ['sigma_a=0.468', 'sigma_b=0.005', 'sigma_c=113.041']\n",
      "MSLE: 0.132, Exp of RMSLE: 1.439, R2 score: 0.881\n",
      "Day of flattening of the infection curve ->  01, May 2020\n",
      "\n",
      "\n",
      "**** Prediction for Japan-JPN as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-02-22 -> 105.0\n",
      "No. of cases on 2020-07-28 -> 29989.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.039', 'sigma_p=0.000']\n",
      "MSLE: 22.404, Exp of RMSLE: 9.670, R2 score: -2.026\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.032', 'c=24381.960'] ['sigma_a=1.342', 'sigma_b=0.003', 'sigma_c=1071.671']\n",
      "MSLE: 1.062, Exp of RMSLE: 2.803, R2 score: 0.900\n",
      "Day of flattening of the infection curve ->  03, Sep 2020\n",
      "\n",
      "\n",
      "**** Prediction for Australia-AUS as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-10 -> 100.0\n",
      "No. of cases on 2020-07-28 -> 14935.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.038', 'sigma_p=0.001']\n",
      "MSLE: 24.520, Exp of RMSLE: 9.670, R2 score: -5.050\n",
      "~~ Logistic Model ~~\n",
      "['a=9.028', 'b=0.012', 'c=32233.665'] ['sigma_a=8.603', 'sigma_b=0.003', 'sigma_c=29684.787']\n",
      "MSLE: 0.623, Exp of RMSLE: 2.202, R2 score: 0.741\n",
      "Day of flattening of the infection curve ->  14, Aug 2021\n",
      "\n",
      "\n",
      "**** Prediction for New Zealand-NZL as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-23 -> 102.0\n",
      "No. of cases on 2020-07-28 -> 1207.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.024', 'sigma_p=0.001']\n",
      "MSLE: 18.770, Exp of RMSLE: 9.670, R2 score: -22.147\n",
      "~~ Logistic Model ~~\n",
      "['a=5.567', 'b=0.213', 'c=1159.461'] ['sigma_a=0.323', 'sigma_b=0.006', 'sigma_c=2.736']\n",
      "MSLE: 0.005, Exp of RMSLE: 1.072, R2 score: 0.985\n",
      "Day of flattening of the infection curve ->  18, Apr 2020\n",
      "\n",
      "\n",
      "**** Prediction for Canada-CAN as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-12 -> 103.0\n",
      "No. of cases on 2020-07-28 -> 114597.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.057', 'sigma_p=0.001']\n",
      "MSLE: 33.354, Exp of RMSLE: 9.670, R2 score: -2.708\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.044', 'c=112813.015'] ['sigma_a=0.805', 'sigma_b=0.002', 'sigma_c=1695.070']\n",
      "MSLE: 1.210, Exp of RMSLE: 3.005, R2 score: 0.973\n",
      "Day of flattening of the infection curve ->  30, Jul 2020\n",
      "\n",
      "\n",
      "**** Prediction for Iran-IRN as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-02-27 -> 139.0\n",
      "No. of cases on 2020-07-28 -> 293606.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.055', 'sigma_p=0.000']\n",
      "MSLE: 38.182, Exp of RMSLE: 9.670, R2 score: -1.986\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.024', 'c=330197.934'] ['sigma_a=0.614', 'sigma_b=0.002', 'sigma_c=20561.105']\n",
      "MSLE: 1.068, Exp of RMSLE: 2.811, R2 score: 0.958\n",
      "Day of flattening of the infection curve ->  10, Nov 2020\n",
      "\n",
      "\n",
      "**** Prediction for Turkey-TUR as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-19 -> 191.0\n",
      "No. of cases on 2020-07-28 -> 227019.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.060', 'sigma_p=0.001']\n",
      "MSLE: 43.861, Exp of RMSLE: 9.670, R2 score: -3.675\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.054', 'c=208673.750'] ['sigma_a=0.946', 'sigma_b=0.002', 'sigma_c=2680.458']\n",
      "MSLE: 0.783, Exp of RMSLE: 2.422, R2 score: 0.966\n",
      "Day of flattening of the infection curve ->  12, Jul 2020\n",
      "\n",
      "\n",
      "**** Prediction for Russia-RUS as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-18 -> 114.0\n",
      "No. of cases on 2020-07-28 -> 818120.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.074', 'sigma_p=0.001']\n",
      "MSLE: 35.339, Exp of RMSLE: 9.670, R2 score: -1.398\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.030', 'c=877346.515'] ['sigma_a=1.172', 'sigma_b=0.003', 'sigma_c=80858.606']\n",
      "MSLE: 4.202, Exp of RMSLE: 7.767, R2 score: 0.908\n",
      "Day of flattening of the infection curve ->  12, Oct 2020\n",
      "\n",
      "\n",
      "**** Prediction for United Kingdom-GBR as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-06 -> 115.0\n",
      "No. of cases on 2020-07-28 -> 300111.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.061', 'sigma_p=0.001']\n",
      "MSLE: 38.971, Exp of RMSLE: 9.670, R2 score: -2.581\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.045', 'c=310983.987'] ['sigma_a=1.136', 'sigma_b=0.003', 'sigma_c=5641.154']\n",
      "MSLE: 2.113, Exp of RMSLE: 4.278, R2 score: 0.953\n",
      "Day of flattening of the infection curve ->  21, Jul 2020\n",
      "\n",
      "\n",
      "**** Prediction for France-FRA as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-01 -> 100.0\n",
      "No. of cases on 2020-07-28 -> 183079.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.056', 'sigma_p=0.001']\n",
      "MSLE: 37.960, Exp of RMSLE: 9.670, R2 score: -3.601\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.055', 'c=166804.289'] ['sigma_a=1.059', 'sigma_b=0.003', 'sigma_c=1840.738']\n",
      "MSLE: 1.355, Exp of RMSLE: 3.203, R2 score: 0.961\n",
      "Day of flattening of the infection curve ->  21, Jun 2020\n",
      "\n",
      "\n",
      "**** Prediction for Italy-ITA as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-02-24 -> 132.0\n",
      "No. of cases on 2020-07-28 -> 246286.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.054', 'sigma_p=0.001']\n",
      "MSLE: 43.168, Exp of RMSLE: 9.670, R2 score: -3.856\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.056', 'c=243030.882'] ['sigma_a=0.895', 'sigma_b=0.002', 'sigma_c=2068.458']\n",
      "MSLE: 1.156, Exp of RMSLE: 2.931, R2 score: 0.973\n",
      "Day of flattening of the infection curve ->  13, Jun 2020\n",
      "\n",
      "\n",
      "**** Prediction for Brazil-BRA as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-15 -> 121.0\n",
      "No. of cases on 2020-07-28 -> 2442375.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.080', 'sigma_p=0.000']\n",
      "MSLE: 32.941, Exp of RMSLE: 9.670, R2 score: -0.713\n",
      "~~ Logistic Model ~~\n",
      "['a=9.998', 'b=0.025', 'c=2159468.826'] ['sigma_a=2.364', 'sigma_b=0.006', 'sigma_c=637699.122']\n",
      "MSLE: 5.820, Exp of RMSLE: 11.161, R2 score: 0.735\n",
      "Day of flattening of the infection curve ->  19, Nov 2020\n",
      "\n",
      "\n",
      "**** Prediction for India-IND as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-17 -> 125.0\n",
      "No. of cases on 2020-07-28 -> 1483156.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.076', 'sigma_p=0.000']\n",
      "MSLE: 27.423, Exp of RMSLE: 9.670, R2 score: -0.560\n",
      "~~ Logistic Model ~~\n",
      "['a=9.994', 'b=0.024', 'c=1066388.221'] ['sigma_a=3.159', 'sigma_b=0.008', 'sigma_c=441540.841']\n",
      "MSLE: 6.032, Exp of RMSLE: 11.657, R2 score: 0.650\n",
      "Day of flattening of the infection curve ->  27, Nov 2020\n",
      "\n",
      "\n",
      "**** Prediction for United States-USA as of 2020-07-29 ****\n",
      "Initial no. of cases on 2020-03-03 -> 103.0\n",
      "No. of cases on 2020-07-28 -> 4290263.0\n",
      "~~ Exponential Model ~~\n",
      "['p=0.079', 'sigma_p=0.001']\n",
      "MSLE: 45.859, Exp of RMSLE: 9.670, R2 score: -1.501\n",
      "~~ Logistic Model ~~\n",
      "['a=10.000', 'b=0.024', 'c=4262579.365'] ['sigma_a=1.038', 'sigma_b=0.003', 'sigma_c=503190.213']\n",
      "MSLE: 5.149, Exp of RMSLE: 9.670, R2 score: 0.903\n",
      "Day of flattening of the infection curve ->  17, Nov 2020\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "iso_codes = ['CHN','THA','KOR','JPN','AUS','NZL','CAN','IRN','TUR','RUS','GBR','FRA','ITA','BRA','IND','USA']\n",
    "\n",
    "for iso_code in iso_codes:\n",
    "    cases_starting_point = 100\n",
    "    \n",
    "    df = data[(data.iso_code == iso_code) & (data.total_cases >= cases_starting_point)][['date','total_cases','new_cases']]\n",
    "    df['time'] = df.reset_index().index\n",
    "\n",
    "    # cutoff_date = '2020-05-31'\n",
    "    # df = df[df.date <= cutoff_date]\n",
    "\n",
    "    x = np.array(list(df.time.values))[0:-1]  \n",
    "    y = np.array(list(df.total_cases.values))[0:-1]\n",
    "\n",
    "    y0 = y[0]\n",
    "    yf = y[-1]\n",
    "\n",
    "    country_name = data[data.iso_code == iso_code].country_name.values[0]\n",
    "\n",
    "    print(\"**** Prediction for {}-{} as of {} ****\".format(country_name, iso_code, df.date.max()))\n",
    "    print(\"Initial no. of cases on {} -> {}\".format(df.date.min(), y0))\n",
    "    print(\"No. of cases on {} -> {}\".format(df.date.values[-2], yf))\n",
    "\n",
    "    # Exponential Model\n",
    "    fit_exp = curve_fit((lambda tt,pp: expo_model(tt,pp,N0=y0)),x,y,p0=[0.10])  ##Take N0 from initial number of reported cases\n",
    "    pi, sigma_pi = fit_exp[0][0], np.sqrt(fit_exp[1]).flatten()[0]\n",
    "    exp_y_pred = expo_model(x,pi)\n",
    "    exp_MSLE=sklm.mean_squared_log_error(y,exp_y_pred)\n",
    "    exp_RMLSE = np.exp(np.sqrt(MSLE))\n",
    "    exp_r2_score = sklm.r2_score(y,exp_y_pred)\n",
    "    print(\"~~ Exponential Model ~~\")\n",
    "    print([ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('p','sigma_p'),(pi,sigma_pi))])\n",
    "    print(\"MSLE: {:.3f}, Exp of RMSLE: {:.3f}, R2 score: {:.3f}\".format(exp_MSLE, exp_RMLSE, exp_r2_score))\n",
    "\n",
    "    # Logistic Model\n",
    "    fit_i = curve_fit(logistic_model,x,y,p0=np.random.exponential(size=3), bounds=([0,0,0],[10,10,1000000000000]))\n",
    "    ai,bi,ci = fit_i[0]\n",
    "    sigma_ai, sigma_bi, sigma_ci = np.sqrt(np.diag(fit_i[1]))\n",
    "    log_y_pred = logistic_model(x,ai,bi,ci)\n",
    "    log_MSLE=sklm.mean_squared_log_error(y,log_y_pred)\n",
    "    log_RMLSE = np.exp(np.sqrt(log_MSLE))\n",
    "    log_r2_score = sklm.r2_score(y,log_y_pred)\n",
    "    print(\"~~ Logistic Model ~~\")\n",
    "    print([ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('a','b','c'),(ai,bi,ci))], [ss+'='+'{:.3f}'.format(xx) for ss, xx in zip(('sigma_a','sigma_b','sigma_c'),(sigma_ai,sigma_bi,sigma_ci))])\n",
    "    print(\"MSLE: {:.3f}, Exp of RMSLE: {:.3f}, R2 score: {:.3f}\".format(log_MSLE, log_RMLSE, log_r2_score))\n",
    "    \n",
    "    perc_flat = 0.98\n",
    "    sol = int(fsolve(lambda x : logistic_model(x,ai,bi,ci) - perc_flat*int(ci), bi))\n",
    "    datesol = datetime.strftime(datetime.strptime(df.date.min(),'%Y-%m-%d') + timedelta(days=sol), ' %d, %b %Y' )\n",
    "    print('Day of flattening of the infection curve -> {}'.format(datesol))\n",
    "    print('\\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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}