{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sat Feb 15 20:36:55 JST 2020\n"
     ]
    }
   ],
   "source": [
    "!date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as ss\n",
    "import scipy.special as sps\n",
    "import arviz as az\n",
    "\n",
    "# Make inline plots raster graphics\n",
    "from IPython.display import set_matplotlib_formats\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from termcolor import colored\n",
    "\n",
    "import seaborn as sns\n",
    "clrs = sns.color_palette(\"Spectral\", 6)\n",
    "def set_plot_style(usetex=False):\n",
    "    sns.set_style('white', {'axes.linewidth': 0.5})\n",
    "    sns.set(style='white', font_scale=1.1,#context='paper',\n",
    "            rc={'xtick.major.size': 6, 'ytick.major.size': 6, 'legend.fontsize': 14,\n",
    "                'text.usetex': usetex, 'font.family': 'serif', 'font.serif': ['Verdana'],\n",
    "                'text.latex.preamble': r\"\\usepackage{type1cm}\"}) \n",
    "    plt.rcParams['xtick.major.size'] = 6\n",
    "    plt.rcParams['xtick.major.width'] = 1\n",
    "    plt.rcParams['ytick.major.size'] = 6\n",
    "    plt.rcParams['ytick.major.width'] = 1\n",
    "    plt.rcParams['xtick.bottom'] = True\n",
    "    plt.rcParams['ytick.left'] = True\n",
    "set_plot_style()\n",
    "    \n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "func_dict = {\"q2.5\": lambda x: np.percentile(x, 2.5), \n",
    "             \"q25\": lambda x: np.percentile(x, 25), \n",
    "             \"median\": lambda x: np.percentile(x, 50), \n",
    "             \"q75\": lambda x: np.percentile(x, 75), \n",
    "             \"q97.5\": lambda x: np.percentile(x, 97.5)}\n",
    "\n",
    "def get_stats(cmdstan_data):\n",
    "    # include mean and hpd\n",
    "    stats = az.summary(cmdstan_data,credible_interval=0.95).loc[:, ['mean','hpd_2.5%','hpd_97.5%','ess_bulk','ess_tail','r_hat']].reset_index().rename(columns={'index':'var', 'hpd_2.5%':'hpd2.5', 'hpd_97.5%':'hpd97.5'})\n",
    "    stats = az.summary(cmdstan_data,credible_interval=0.50).loc[:, ['hpd_25%','hpd_75%']].reset_index().rename(columns={'index':'var', 'hpd_25%':'hpd25', 'hpd_75%':'hpd75'}).\\\n",
    "        merge(stats, left_on='var', right_on='var')\n",
    "    # include percentiles\n",
    "    stats = az.summary(cmdstan_data, stat_funcs=func_dict, extend=False).reset_index().rename(columns={'index': 'var'}).merge(stats, left_on='var', right_on='var')\n",
    "    stats['time'] = stats['var'].apply(lambda st: st[st.find(\"[\")+1:st.find(\"]\")])\n",
    "    stats['time'] = ['NA' if \"[\" not in y else int(x)+1 for x,y in zip(stats['time'],stats['var'])]\n",
    "    stats['var'] = stats['var'].apply(lambda st: st[:st.find(\"[\")] if \"[\" in st else st)\n",
    "    return stats.loc[:,['var','time','mean','hpd2.5','hpd25','hpd75','hpd97.5','q2.5','q25','median','q75','q97.5','ess_bulk','ess_tail','r_hat']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "resultsdir = \"../../../Hokkaido_Backup/Wuhan Serial Interval 2020/stan-sims-certain_and_probable\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "import scipy.special as sps\n",
    "\n",
    "def find_params_Weibull(p,mu,SD):\n",
    "    θ, k = np.exp(p)\n",
    "    return (θ*sps.gamma(1+1/k)-mu, θ*np.sqrt(sps.gamma(1+2/k)-sps.gamma(1+1/k)**2)-SD)\n",
    "\n",
    "def weibull_pdf(x,θ,k):\n",
    "    return (k / θ) * (x / θ)**(k - 1) * np.exp(-(x / θ)**k)\n",
    "\n",
    "weibull_cdf = lambda x,θ,k: -np.expm1(-(x/θ)**k)\n",
    "\n",
    "mean_SARS = 8.4\n",
    "sd_SARS = 3.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[31mlognormal-truncated\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "folder = \"lognormal-truncated\"\n",
    "distrib, truncation_type = folder.split(\"-\")\n",
    "\n",
    "print(colored(folder, 'red'))\n",
    "posterior_glob = !cd \"{resultsdir}/{folder}\"; ls trace-*\n",
    "cmdstan_data = az.from_cmdstan(posterior = [resultsdir+\"/\"+folder+\"/\"+x for x in posterior_glob], \n",
    "                               log_likelihood=\"log_likelihood\")\n",
    "param1 = cmdstan_data.posterior.param1.values.ravel()\n",
    "param2 = cmdstan_data.posterior.param2.values.ravel()\n",
    "\n",
    "xstep = .01\n",
    "xmax = 20\n",
    "ymax = 0.12\n",
    "x = np.arange(.01,xmax+xstep,xstep)\n",
    "\n",
    "nsamples = 10000\n",
    "# a bit faster than to use scipy.stats functions\n",
    "lognormpdf = lambda x, mu, sigma: (np.exp(-(np.log(x) - mu)**2 / (2 * sigma**2)) / (x * sigma * np.sqrt(2 * np.pi)))\n",
    "lognormcdf = lambda x, mu, sigma: 0.5 + 0.5*sps.erf((np.log(x) - mu)/np.sqrt(2)/sigma)\n",
    "yy = np.stack([[lognormcdf(xx+xstep,param1[idx],param2[idx])-lognormcdf(xx,param1[idx],param2[idx]) for xx in x] for idx in range(nsamples)]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAD4CAYAAABsdWSLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd1gU1/v271k6SAfpRVCKFEVABAEVxIbYEDVRgzXRWKNR7L0nxtgjfO2KGo2KgKIURUFsEQuiIkpXBJG2dJZ9//Dd/UnYMouwtPO5Lq7IzJmZZwh7c855GsVms9kgEAgEgkgwWtoAAoFAaIsQ8SQQCIRGQMSTQCAQGgERTwKBQGgERDwJBAKhERDxJBAIhEYgVvHMy8vDjBkzYGtrCzc3N5w6darRY6OiouDl5QUbGxv4+fkhJSWFe+769esYM2YM7O3t4eDggOnTpyMjI6PZ3otAIHQ8xCqeq1atgq6uLuLj43HgwAHs27cPiYmJIo/Nzc3F0qVLsWrVKty7dw9eXl5YsGABOCGrCgoKWL16NeLi4nD79m1oaWnh999/F9t7EgiE9o/YxLO0tBRxcXFYtGgRFBUVYWNjA29vb1y7dk3ksVFRUXByckK/fv2goKCA6dOn49OnT3j9+jUAwNXVFXZ2dpCTk0NZWRmKiopgY2MjrlclEAgdALGJZ0ZGBuTl5aGiosI9ZmxsjPT0dJHHpqWlQV9fn3tOQkICBgYGSEtLq3efyZMnw9XVFcrKypgxYwZtW2tra5GdnY3a2lra1xAIhI6F2MSzsrISsrKy9Y7JysqivLxc5LEVFRWQk5MTeq+TJ08iISEBFRUVWL9+PU+79u7dC3Nz83pfVlZW8PT0RG5ursjvSSAQOgZiE09ZWdkGM7nKysoGIkhnrJycHGpqamjdS01NDdOnT0dERARPu+bNm4fXr1/X+4qOjhbp3QgEQsdDbOJpZGSEkpISfPr0iXssLS0NRkZGIo81NjZGamoq9xyLxUJmZiaMjY15PrusrAzS0tJN9CYEAoEgRvFUVFSEs7Mz9uzZg9LSUjx//hxhYWEYMmQIAMDf358bjiRsrIeHBx4+fIjY2FiUlZUhMDAQ6urqsLS0BABs2LABd+/eRVlZGTIzM/Hnn39i+PDh4npVAoHQAZAU58M2bdqEFStWwMXFBcrKypgzZw4cHBwAAFlZWSgsLKQ1Vk9PD1u3bsXGjRuRm5sLCwsL7N69GxRFAQDU1dWxfv165ObmQlVVFcOHD8e8efPE+aoEAqGdQ5F6ng3Jzs6Gp6cnoqOj63n1CQQCgQNJzyQQCIRGQMSTQCAQGgERTwKBQGgERDwJBAKhERDxJBAIhEZAxJNAIBAaARFPAoFAaAREPAkEAqEREPEkENow+fn52Lx5MwYOHAhra2u4uLhg4sSJOHXqFKqrqwEA9+7dw4wZM+Do6AgbGxv4+PggMDCQW1yHyWSiR48eOHHiBM9nbN++nZve7OHhUa8CWXZ2NgBg2bJl3GPW1tZwdXXFzJkzERYWJtL7sNlsvH37FoMHD8bhw4cbnK+oqMAff/wBT09P2NnZYcKECUhISBDpGU2FWNMzCQRC05GdnY0JEyZAT08PS5YsgampKZhMJuLi4rBv3z7Y29vjzZs3WLZsGSZOnIj58+dDQUEBjx8/xp9//okHDx4gMDAQnTp1wsCBAxESEoIffvih3jPq6uoQHh7OPR4cHIzDhw/j9u3bOHz4MLS1tblj+/Tpg82bN6O2thYfPnxAZGQkli9fjoSEBGzevJnWO3Xv3h11dXV8z8+bNw8FBQXYtGkT1NXVcePGDcyYMQPHjx/npm+LDTahAVlZWWwzMzN2VlZWS5tCIPDlp59+Yo8aNYpdVVXV4FxBQQE7JyeH7eDgwN68eXOD88+ePWObm5uzL126xGaz2ew7d+6wzczM2KmpqfXGxcfHsy0tLdm5ubncY3v27GF7e3vXGxcQEMD+8ccfGzwnMTGRbWVlxX2OMFJTU9mpqalsV1dX9v/+97965969e8c2MzNjP336tN7x+fPns2fNmkXr/k0JWbYTCG2Q4uJixMbGYtq0aTzLLaqpqSEpKQklJSU8uyjY2NjA2dmZW+fWxcUFWlpauHLlSr1xISEh3HONoWfPnhg+fDjOnj1La7ypqSlMTU0hJSXV4FxxcTGAL1XXvsbJyQnJycmNsu9bIMt2AuE/TJ48GTExMWJ9poeHB06ePEl7fEZGBurq6tC9e3e+Y9LT09GpUyd07tyZ53lTU1PcvXsXAMBgMDBixAhcuXIFCxcuBEVRqKiowI0bN7Bp0ybRXuY/WFtbIzIy8pvuAQBmZmZQUlLC/v37sXz5cqipqeHjx49ITk7m2ZGiuSEzTwKhDcKmWQyNU6aRzvnRo0fj/fv3ePDgAYAvjRYZDAY8PT0bb+h/npOTkwMbG5t6Xzk5ObTuIS8vjwMHDuDly5dwcXFB9+7dMWbMGKSkpDSYjYoDMvNsxTCZTPz777/4/PkzampqoKmpCXd3d0hISLS0ae0aUWaALYWBgQEoikJSUhJMTU15jjEyMkJpaSny8vJ4zj7fvn1br5ODqakpbG1tERISAicnJ4SEhGDo0KEN+omJSnJyMszNzQEAnTt3xuXLl+ud5zcz5oWjoyPCw8ORn58PFouFzp07Y8mSJY3eVvgWyMyzFZKXl4e///4bf//9NwBAV1cXhoaGKC0txY4dO5CZmdnCFhJaGjU1NTg5OWH//v2orKxscL6kpAQuLi5QVFREYGBgg/PPnj1DQkICtzsDh1GjRuH69evIzs7G3bt3MWrUqG+yMykpCaGhoZg4cSIAQEpKiruvKWh/UxiamprQ1tZGeno6bty4gZEjR36TnY2BzDxbGQ8fPkRGRgb09fVBURQUFRW5M01lZWXo6enh/Pnz8PHxgZmZWQtbS2hJ1qxZg++//x6+vr6YM2cOzMzMUFFRgfv37+P48eM4c+YMVq1aheXLl6OiogIjR46EoqIiN1TJ3d29QXsab29vbN26Fb/++iv09PQahP/k5uaipKQENTU1yM7Ohra2NiQlv8hIZWUlsrOzUVlZifz8fERHR+PChQsYO3Yshg0bRuudSkpKAHwJkaqqqkJJSQlkZGQgIyMDAPjw4QPy8/MhLS2NZ8+eYffu3RgyZAgGDhz4rT9OkSGV5HnQUpXki4uLERERAQMDAygpKfFdntfU1CA8PByzZ8+u19ue0PHIycnBwYMHERcXh/z8fMjJycHa2ho+Pj4YOXIkJCUlkZCQgMDAQDx79gxVVVUwNDTEqFGjMHXqVJ6zvvnz5+P69euYN28e5s6dW++ch4dHvT1Kzmdk2bJluHTpEgBAUlISysrKsLW1xbhx4+Dh4UH7fTjL+6+ZO3cut41OQkICZs6cCSkpKZiYmMDPzw/jxo0DgyH+RTQRTx60hHjW1dXh1KlTMDc3h4KCgtB9zeLiYsTFxWHu3LnfvCdFIBBEh+x5thLCwsJgaGgIKSkpWg4hZWVlWFpaIiQkRGBGBoFAaB6IeLYCkpOTISEhAQUFBcjJydG+zsTEBO/fv0deXl4zWkcgEHhBxLMVkJSUBHV1dSgoKIh8rbOzM8LCwmjH/REIhKaBiGcLk5WVBUlJScjIyAgNaOaFsrIyKioqkJWV1QzWEQgEfhDxbGHi4+Ohp6fHDcVoDH379sWVK1e4JcYIBELzQ8SzBXnx4gVkZWUhKSnZqFknBzk5ObBYLLL3SSCIESKeLcirV6+gpaUlkpOIH05OTrh27RpYLFYTWEYgEIRBxLOFKCoqAovFgrS09DfNOjkoKyuDyWRyy3YRCITmhYhnCxEXFwdDQ8MmmXVy6NatG+Li4ojnnUAQA0Q8W4jS0lIwGIwmmXVyMDExwbt371BWVtZk9yQQCLyhJZ6cRlKEpuHNmzdQVlbmWQH8W+nUqRNSUlKa/L6E1kdtbS1OnDiBESNGwM7ODk5OThgxYgQ2bdrEbczGYcyYMTA3N29QDg4A9u7dy23eZmlpiT59+mDp0qXIzc2tN668vBz79+/H0KFD0aNHDzg7O2PMmDH4/fffUVhY2Kzv2hqhVVXJ1dUVI0aMwJgxYwRWribQ4+XLl9DW1v6m8CR+ODk54ebNm+jevTvJeW/nbNiwAZGRkQgICICNjQ2qq6uRlJSEixcv4unTp9y6DG/evEFycjKcnZ1x6dIlnmXmjI2NcfjwYbBYLKSnp2PLli2YNWtWPbGdN28esrKysHjxYnTr1g3l5eVITEzE33//jYEDB0JVVVVs794aoCWe27dvx+XLlzFhwgR06dIFvr6+8PHx6XA/rKaCyWQ2270lJSXBZrNRWFgIHR2dZnsOoeXhtMz4WgwtLS3h5+dXb+vm0qVLcHR0xKxZszBlyhS8f/8eurq69e4lJSXFFVsjIyNUVFRgwYIF+PTpEzQ0NPD582fExcVh//799cq/WVtb4/vvv++Qq1Nay/YBAwZg9+7diIuLw7hx4xAaGgp3d3fMnz8fsbGxpDCFCGRkZEBGRqZJHUX/xcLCAgkJCcRx1M7R0tLCvXv3UFFR0eAcJ9WXxWLhypUrGDlyJJycnKCjo4OQkBCh966qqoKEhAS3vYWioiLk5eVx586dBskYEhISzfr73FoRqRiykpISJk6cCHd3d/zxxx+4du0abty4AU1NTYwaNQq+vr4wNjZuJlPbB5zlFKeAbHNgaGiIyMhIlJWVoVOnTs32nPaMnp4e33Pbt2/HpEmTAACnTp1CQEAA37Ff174cMmQInj9/LnQcXdasWYOFCxeiT58+sLa2hoWFBWxtbeHp6cn9/x4XF4fS0lIMGTIEFEVhxIgRuHz5MmbPns3znnV1dXj16hUOHjwIX19f7taSlJQU1q5dizVr1iAsLAw2NjawsLCAnZ0d+vfv3yxbUK0d2t728vJyXLx4EZMnT8bgwYORnZ2NNWvW4O7du1i5ciVevHhBu1p0R6aqqqpJPez8oCgKBQUFzf4cQsvRt29f3Lx5E9u2bUPPnj2RkZGBDRs2wMvLCy9fvgQAXL58GQMHDuSK6ciRI5Geno4nT57Uu1dqaipsbGxgZWWF0aNHw93dHWvWrKk3ZtSoUbh16xZWr14NMzMzvHz5EkuWLIG3tzc+fPggnpduRdAqhrx06VJERkZCQUGB6zjq2rVrg3FZWVkwMDBoFkPFSXMVQ66pqcG5c+fQvXv3Zl/mfPjwAXl5efDz8yMN4zoQTCYT48ePh4WFBdatW4e+ffuiqqqqwbgJEyZg/fr1AL5428PCwnDgwAFUVVVh/vz5MDIywl9//SW0v1B+fj5GjBgBPz8/LFq0qFneqbVCa+1YWVmJXbt2wd3dXWC5+/YgnM3JixcvoKGhIZYljo6ODpKSklBaWkpadbRDqqur8fHjxwafuU6dOkFPTw/V1dW4evUqFBQUcO7cuXpjIiIicObMGaxcuZIbLsdpzAYAQUFBmDBhAlauXInt27eDoigUFRWhrq4Oampq9e6lqakJNTW1DukwoiWeZmZmcHBwaCCcqampePjwIb777rtmMa69kZmZCW1tbbH1W5GQkEB+fj4Rz3ZIZWUlfHx8MGrUKHh4eMDAwAClpaWIiIhAQkICjh8/jh07dmDo0KGwtLSsd62GhgaCgoIQHR2NoUOHNrh3ly5dcPDgQUyZMgWdO3fGr7/+ig8fPsDf3x9jx46Fu7s7tLW18fnzZ5w/fx55eXkdUgNofYr379/P7Wr3NTU1Nfjtt9+a3Kj2Cq8Wsc2JtbU1Hjx4gNraWrE+l9D8KCgoICAgABkZGVi9ejV8fHwwbdo0vHv3DqdPn4a6ujoSExPh4+PT4FpNTU306dOHZ8A8h169euH333/H4cOHcfLkSRgaGmLWrFl48uQJfv31V3h7e2P27Nmorq7GuXPn6vV/7yjQ2vO0sLDAzZs368UN1tbW4ujRozh+/Dji4uKa1Uhx01x7nkeOHIGdnZ1Yg9cjIyPxww8/kNkngdDECFy2W1hYgKIoUBTFs32ohIQEli9f3mzGtScKCgogKSnZLCmZgpCQkEBBQQERTwKhiREonidOnACbzYa/vz/++OMPaGhocM/JyMjA0NCQZBnRJDExEV26dBF7f2lbW1vcvXsXxsbGxOtOIDQhAsWzd+/eAL4U7SV8G6WlpVBSUhL7czU0NPD48WMwmUwoKyuL/fkEQntFoHju27cPU6ZMQVRUlMCb8Co0QKhPZWWlWILjeSEpKYn8/HwingRCEyJQPC9evAg/Pz/s2bOH7xiKooh4CqGurg5lZWVi3+/kYGdnh/j4eJiYmIh924BAaK8IFM+YmJh6/yU0jjdv3qBz585CszWaC1VVVRQVFYHJZLbI1gGB0B6hNQ0pKCjAixcvuN+npqbi6NGjuH//frMZ1p5ISUmBtrZ2iy3bgS9Ld5LrTiA0HbTEc+PGjTh16hQAIDc3F35+fjh79izmzJmD8+fPN6uB7YHW0E/d1taWlKkjEJoQWuL5+PFjDB8+HMCXvFgDAwNcv34dO3bswJEjR5rVwPZARUVFi+81amhooKCggGftRwKBIDq0PtFlZWXcytOJiYlwcnIC8KXhWEcsRSUKTCYTbDa7xfY7v4bNZqO0tLSlzSAQ2gW0xNPKygoXLlzA48ePcefOHbi6ugIA3r59C01NzWY1sK3z7Nkz6OvrtwrxNDExwb///tvSZhAI7QJa4rl8+XKEhYXh+++/h4ODA9zd3QEAgYGBXCEl8KagoAAKCgrf5CxisVgoLy//5v1KU1NTvH37VuwFSgiE9gitknSWlpa4desWCgsLufX86urq8McffzSo70eoT21tbaNFr7KyEsHBwYiIiMDnz5+hrKwMPz8/jBo16puKi5SUlJDOmgTCN0Lbi0FRVD2hZDAY0NPT65CNn0Shsc6iDx8+YO7cuQgODsbnz58hJSWF4uJi/O9//8PPP/+MrKysRtmjqalJ+roTCE0ArU/127dvMW3aNDg6OsLS0rLBF13y8vIwY8YM2Nraws3NjRv+1JixUVFR8PLygo2NDfz8/OoJQkxMDMaOHQt7e3v06dMHK1asQHl5OW07m4qKigrU1taK3OytpKQEy5cvR3p6OvT09LBz505cvXoV27Ztg6GhITIzM7FgwQKkpaWJbJO1tTWePHlCanwSCN8IrU/1kiVLoKysjI0bN0JVVbXR+3erVq2Crq4u4uPjkZ6ejpkzZ8LKygp2dnYijc3NzcXSpUuxa9cuODg44PTp01iwYAGuXr0KiqJQVlaG2bNnw8nJCZWVlfjll1+wf/9+LFmypFF2N5bXr19DS0tLJGcRm83Gzp07kZ2dDRMTE+zatYvbRtbBwQH79+/Hxo0b8eDBA6xZswYHDhzgtoelg4SEBNhsNkpKSsiWC4HwDdASz8zMTAQHB8PMzKzRDyotLUVcXBzu3r0LRUVF2NjYwNvbG9euXWsgnsLGRkVFwcnJCf369QMATJ8+HUFBQXj9+jUsLCzqVc/u1KkThgwZ0iIppjk5OdDQ0BBp2X737l3Ex8dDXl4eGzdu5AonBzk5Oaxduxa//PILUlJSsHPnTqxdu1akP2gKCgrIysoi4kkgfAO0PtU2NjbIzMz8pgdlZGRAXl6+XlFeY2NjpKenizw2LS2tXoV3CQkJGBgY8F3GPn36FCYmJjzP7d27F+bm5vW+PD09G/GGDamtrRVJ1FgsFoKCggAAU6dOhZaWFs9xMjIyWLVqFeTl5REXFyfyHwY7OzskJCSAxWKJdB2BQPg/aM08+/btix07dvAtaebo6Cj0HpWVlQ08vLKysjz3IoWNraioqFeYWdC9YmJiEBsbi0uXLvG0a968eZg3b169Y5w2HN+KqGmZd+/eRXZ2NrS1tXn2nvkaXV1dzJ49Gzt37sTBgwfh4OBAu+ScrKwsqqurUVZWRgqFEAiNhJZ4/v777wCAyZMnNzhHURRevnwp9B6ysrINnBSVlZU8vfXCxsrJyTUQJl73unfvHpYvX469e/dyM6TESWVlpUjV269cuQIA8PX1peVkGjJkCKKiovD06VMcPnxYpL7Z0tLS+PjxIxFPAqGR0Fq2v3r1iu8XHeEEACMjI5SUlODTp0/cY2lpaTy77gkba2xsjNTUVO45FouFzMxMGBsbc49FRUVh0aJF2L9/P7civjipqKgAi8Wi7SzKzMxEYmIiZGVl4eXlResaiqKwYMECMBgMRERE8NwC4UevXr1w+/ZtUiiEQGgkIgUgstlsfPz4sVF7ZYqKinB2dsaePXtQWlqK58+fIywsDEOGDAEA+Pv7c8ORhI318PDAw4cPERsbi7KyMgQGBkJdXZ0bNnX27Fls2bIFR44cgYODg8i2NgVpaWnQ0NCgHaZ0/fp1AMCAAQPQqVMn2s8xNDTE8OHDUVdXh8DAQNrXKSoqoqysDGVlZbSvIRAI/wct8SwqKsIvv/wCGxsbDBgwgOs8mjNnDvbt20f7YZs2bUJWVhZcXFwwe/ZszJkzhytuWVlZKCwspDVWT08PW7duxcaNG+Hk5ITo6Gjs3r2b65wJDw/Hhw8fMGbMGHTv3p379eDBA9q2fisZGRlQU1Oj5TBis9m4desWANCedX7N5MmTIS8vjwcPHuD58+e0r5OQkEBeXp7IzyMQCDT7ti9evBj5+flYvHgxJk2ahCtXrqBLly6IiorCH3/8gatXr4rDVrHRFH3bQ0JCoKOjQysG88WLF1iwYAE0NTVx+vTpRmUkHTt2DKdOnYKjoyO2bt1K65rPnz/j3bt3mDRpUosWaiYQ2iK0PqVxcXEICAhAjx496n3ITE1NkZOT02zGtWVEyeCJjY0FAPTv37/RdT9Hjx4NWVlZPHz4EK9fv6Z1jZqaGgoLC1sk+4pAaOvQ/qTW1dU1OJaeni5SdktHoqqqivZsjtPOpG/fvo1+nrKyMkaMGAEACA4Opn0dRVH1tksIBAI9aInn4MGD8eeff4LJZAL48oF79eoVtm/fjsGDBzergW0RNpuNqqoqWp727Oxs5OTkQFFRUaQ6AbwYO3YspKSkEB8fj3fv3tG6xsrKCnfv3v2m5xIIHRFa4rlixQqoq6ujb9++qK6uhq+vL0aPHg1zc3MsXry4uW1sc+Tl5UFeXp6Wp/3hw4cAvuStixITygs1NTUMGzYMwJeIAzro6OggNzeX1PgkEESEVhyNrKwsduzYgfnz5yM1NRUsFgtdu3blGaNJ+NJdtHPnzrT2LzlL9qaKRR03bhxCQ0MRGxuLH3/8sUEmFj9IjU8CQTRE8k7o6+ujf//+8PT0JMIpgKKiIlrV4ysqKvD06VNQFNVk8ahaWlpwdXUFi8XiZiwJo2vXrmIN4yIQ2gN8Z56TJ0+m7fA4ceJEkxnUHqDraX/58iVqampgZmYGVVXVJnv+mDFjcPv2bYSFhWHixImQkZERON7Y2BhRUVGoqqoSOpZAIHyB78zTyckJvXv3Ru/evcFkMiElJcX9nvP15s0b6OnpidPeNgGvyARevHr1CgDQvXv3Jn2+lZUVzMzMUFJSgujoaKHjOX8kS0pKmtQOAqE9w3fmOXfuXO6/Q0NDsXbtWvTo0aPeGA0NDeKp5QHdakqcugDf6mX/LxRFwdfXF1u3bsXFixcxdOhQoasIQ0NDPHnypFEZTgRCR4TWnueHDx94Lud69erFdXgQ/o+qqiqhnnY2m80VTwsLiya3wd3dHerq6khPT8fjx4+Fju/WrRtevXqF6urqJreFQGiP0BJPOzs7HDt2rEFBkJiYGMjLyzeLYW2V6upq1NbWCg07+vjxI4qKiqCkpNQs5fKkpKS4QfMXL14UOp7BYIDNZqO4uLjJbSEQ2iO0QpU2bdoEf39/eHp6wsbGBpKSkkhNTUVGRgZ27NjR3Da2Kd6/fw8VFRWhM8+vZ53NlVfu7e2N06dP4/79+8jOzhaap6+jo4MXL16gf//+zWIPgdCeoDXzNDAwQEREBBYsWAADAwOoq6vD19cXUVFR3DJxhC9kZGRAXV1daIwnx1nUHEt2DioqKtyK+Pwq6X9N9+7d8ezZM1RVVTWbTQRCe4F2T1xpaWmMHj26OW1pFzCZTGhqagodJw7xBL4UDLl27RquX7+OqVOnCqwVylm6l5SU0HoHAqEj07gSPgS+0Gn6VlNTw+0z39ziaWJiAjs7O1RWVnILLgtCV1dXpJqgBEJHhYhnE0Onyv67d+9QU1MDfX19sfQQGjVqFIAvPZKExaB2794dz58/J0t3AkEIRDybmNraWqF9gcS1ZOfQp08faGpqIicnB4mJiQLHMhgMUBRFvO4EghCIeDYxlZWVtD3tTR0czw8JCQl4e3sDAK18dz09PbJ0JxCEQEs8Dx8+jIKCgua2pc3DqeMpTDw5+53m5ubiMAsAMGzYMEhKSiIhIUFo3yILCwskJSWRgHkCQQC0xPPq1avo168fZs2ahaioqEZ1z+wIFBcXQ1ZWVmCYUlVVFbKzs8FgMGBiYiI229TU1ODm5oa6ujqEh4cLHMuxv6ioSBymEQhtElri+c8//+DSpUswMTHB+vXr4ebmhm3btnFnUIQv5OTkQFVVVWB2UWZmJurq6qCvrw9paWkxWgduxtHVq1eF5t9369YNCQkJ4jCLQGiT0N7z7NatG5YuXYrY2Fhs3boVaWlpGDlyJMaOHYszZ86gtLS0Oe1sE7x//x5qamoCx3DaY3Tp0kUcJtXD2toaxsbGKCwsRFxcnMCxXbp0QXp6OmkORyDwQWSH0aNHjxAREYEHDx5AQ0MDNjY2CA4Ohqura4dvyUHHWZSWlgagZcSToiju7DM0NFToeAaDgc+fPze3WQRCm4SWeGZlZWHv3r3w9PTEjBkzUF5ejj///BOxsbFYu3YtQkNDceTIEVoNz9ozdIogc2ae4tzv/JqBAwdCXl4ez549E9okzs7ODjdv3hQaekUgdERopWd6eXnB2toa06ZNw/Dhw6GsrNxgjL29Pezt7ZvcwLYEHUcaZ+Ypinh++vQJjx49Ql1dHaSkpNCpUyf06dOnUQVF5OXl4eXlhZCQEISGhmLBggV8x6qrq+Pff/9FWVmZwLROAqEjQks8Q0ND0a1bt+a2pc1TU1Mj0NNeWFiIwsJCyMvLQ4OYAhIAACAASURBVEtLi9Y9nz17htLSUnh4eEBfXx9SUlJ49eoVrly5gkGDBkFOTk5kO318fBASEoKoqCjMmDEDCgoKfMfKyckhJydHrGFVBEJbgNayfcSIEXj//n2D43FxcfDw8Ghyo9oq1dXVAj3tX+930pk1fv78GUVFRRg8eDBsbGygrq4OJSUl9O7dG/7+/oiJiWlUGqWxsTFsbW1RUVGBqKgogWMdHBxw8+ZNEp5GIPwHWuLJb89LRUWFOBT+P5wAebriSYf79++jT58+0NXVbTCj1dHRwbhx43Dz5k3aDee+huM4unLlisA9TTk5OVRVVZFoCgLhPwhctu/btw/AFy/t8ePHoaioyD1XU1ODyMhI2NjYNK+FbYTS0lJIS0sLjfEEQKttc1paGgwNDdGtWze+9zQ0NMSgQYNw69Yt9OvXT2j1+q9xdXWFuro6MjIy8PTpU/Ts2ZPvWA0NDbx8+RLOzs60708gtHcEiienPxGbzcaTJ0/qBXXLyMjA1dUV06dPb14L2wjZ2dlQVVUVuOfJEU9DQ0Oh93v79i2GDBkitM2JlZUVSktL8fjxYzg4ONB2IklKSsLb2xsnTpxASEiIQPHs0aMHoqOjYW9vL/bAfgKhtSJQPE+ePAngywzU39+/3syTUJ/s7GxoaGgIHJOVlQXgS2V+QZSXl0NGRoa2U8nJyQlv3rxBXl4e7WuA/2vTER8fj/z8fL4FkDkz2s+fP0NbW5v2/QmE9gytPc+5c+cS4RRCRUWFwFlZcXExioqKICsrK7RK+/Pnz2FtbU07PIiiKIwfPx6JiYkiOZDU1dW5+e5hYWECx9rY2CA2NpbEfBII/x+BM09LS0tER0dj0qRJApeD0dHRTW5YW0OY04Yz6zQ0NBS6tGYymdDW1hbaB+lrpKWlMWzYMMTFxaFPnz60rx05ciRu3bqF8PBwTJw4ke8fAB0dHTx//hxMJpP8ISUQIEQ8t2zZAhUVFcybN09c9rRZWCyWwFkZ3SU7k8mEjIwM1NXVRbbB3Nwc9+/fR2FhIe3rra2tYWJignfv3uHOnTvchnG8UFFRwatXr+Do6CiybQRCe0OgeHIavpHGb8IRFiDPcRYJE8/k5GR079690Rk9vr6+OHLkCPr160crXZaiKIwcORK7du1CSEiIQPF0cHBAVFQUbG1tISMj0yj7CIT2Al/xfPjwIe2bkJnIF/EUFCr09bJdEEwmE1paWiKFHX2NgoICunbtig8fPtDy6gOAh4cHAgMDkZycjDdv3vDNJmMwGGAwGMjLyxP6R4BAaO/wFc/JkyfTugFFUdy2Eh2Z2tpagYLHydDS19cXeJ+6ujqoqKh8ky2DBg3C/v37oa2tTSu0SE5ODkOGDME///yDkJAQ/Prrr3zHOjk54fr165g2bZpIe7IEQnuDr3hympQRhMNkMiEpKclXPFksFj58+ADgi+OFHx8/foSGhobAXHM6SEhIwMHBAampqejevTuta3x8fPDPP/8gJiYGM2fO5Fn8BQAUFRXBZDLBZDLF0vmTQGitkKlDE5CTkwMVFRW+M7H8/HzU1NRAXV1dYCGPlJQUmJmZQVZW9pttcnZ2xsePH1FWVkZrvL6+PhwdHVFdXS00bKlbt25CiykTCO0dgQ6j5cuXIyAgAKdOnRJ4k7lz5zapUW0NYQHyOTk5AABdXV2B96mqqoKWllajSs39F4qiMHz4cERFRaF379609lD9/Pzw8OFDXLx4Eb6+vnxFvEuXLoiMjERVVRVxHBE6LALFMzs7GywWi5umyYum+KC3dcrKygTGPnL2OwWJZ11dHQA0aQyloaEhpKWlUVhYKDT7CfhS/Njc3ByvX7/GtWvX+EZZUBQFZWVlpKamwsrKqsnsJRDaErTSMzn/JfCGI3z84Iinnp4e3zE5OTnQ1dUVmssuKqNHj8bhw4fRv39/oaFLFEXhu+++w7p163D+/Hn4+PjwbStib2+PmJgYWFhYNDoygEBoy4i05/nixQtERkbixo0b3KUoQXh2EZ1le2ZmJgwNDZu8lYmCggIsLS25oVLCcHFxgaGhIfLy8hATE8N3HEdUhfWAJxDaK7TEMzk5GcOGDYOfnx/WrVuHlStXwsvLC4sXLwaTyWxuG1s9tbW1Arcv6Mw8a2tr+Xq4v5WBAwciLS0N1dXVQscyGAxMmDABAHD27FmBs+revXvj2rVrJN+d0CGhJZ7Lli2DpaUl4uPjER8fj4cPHyI8PBzZ2dlYu3Ztc9vY6qmtreXraa+rq+OKp6Awpdra2mbLGacoCvb29nj79i2t8R4eHujcuTMyMzNx9+5dvuOUlJRQXFxMCiUTOiS0xDM9PR0zZ86Eqqoq91iXLl0QEBAgcGnXURDUfqOgoADV1dVQUVHhm3LJZrPBZrObfL/za/r06YPc3FxUVFQIHSspKQk/Pz8AwJkzZwTOLM3MzBAbG9tkdhIIbQVa4mllZcVzz0xBQaHDF8ctLy+HhIQEX/Gks99ZUFAAdXX1Zg37oSgKAwYMQFJSEq1l9tChQ6GiooLXr18jMTGR77guXbogNTWVdjwpgdBeoJXbPmDAAOzduxeqqqr1PnjR0dHo2rVr81rYynn//j2UlZX5iiedMKV3797BzMys2WMmLS0tcefOHVpl5WRlZTF69GgcPXoUZ8+eRa9evXiOoygKOjo6+Pfff+Hu7t4cZhMIrRKRctsnTZrU4FhHj/PMzs6Gmpoa358DZ+YpyFlUUVEBVVVVsfwsR44ciZCQEPTt21dobvqIESNw9uxZPH78GK9fv+bbftjW1hbXrl2Dg4NDs249EAitCZLb/o0wmUx07tyZ73k6M08WiyW2PHEtLS3Iy8ujqKgIampqAscqKirCx8cHf//9N86cOYN169bxHEdRFAwNDZGQkCCwpB2B0J4gue3fiLB+5nTClNhstsCc96Zm1KhRePz4Ma2Wxb6+vpCSkkJ8fDy3JikvrKys8OTJE5SXlzelqQRCq0VghtHXXLp0CYmJiaipqWlwbuvWrU1qVFtCkHiy2WyhDqPCwkIoKSk1STEQunTq1AmGhob4+PGjQFEHvvQ5Gjx4MMLCwnDu3DksWbKE5ziKotClSxfcuXMHgwcPbg6zCYRWBa2Z565du7BhwwYUFRUhNDSUG9cXHR2N4uLiZjWwtSMoiLywsBCVlZVQVFTkuyx/8+YNTE1NxSqewJfOmcnJybQC58eNGwcGg4GoqCh8/PiR7zgLCwskJSWR2SehQ0BLPENCQvDbb79hz549UFRUxM8//4ytW7di9uzZQjtBfk1eXh5mzJgBW1tbuLm5CazWJGxsVFQUvLy8YGNjAz8/P6SkpNQ7X11djStXrjS7B7impkaos0jQfmdFRQXU1dXF7niTlJREz549kZ6eLnSsrq4u+vXrBxaLhQsXLvAdR1EUzMzMSOwvoUNASzwLCgq4IUkqKir49OkTgC81IyMjI2k/bNWqVdDV1UV8fDwOHDiAffv28Y0hFDQ2NzcXS5cuxapVq3Dv3j14eXlhwYIF3DCqwsJC2NnZISAggLZtjUVQBfnW5iz6L25ubsjOzkZlZaXQsd999x0A4OrVqygqKuI7ztTUFCkpKY2efVZWVqKmpoakfBJaPbTEU0VFhSsEPXv25Armp0+faDkdAKC0tBRxcXFYtGgRFBUVYWNjA29vb1y7dk3ksVFRUXByckK/fv2goKCA6dOn49OnT3j9+jUAQFVVFS9evMCxY8do2dZY2Gy2wMZvdMKUALRYeA9FUXB3d8eLFy+EipWJiQmcnJxQVVWFixcvCrynpaUlbty4Ue94fn4+oqKisHv3bixevBhjx46Fs7Mzhg0bVm9c7969YWxsjK5du8LJyQne3t6YOXMmtm3bhkePHjX+ZQmEJoaWw8jNzQ0vXryAi4sLfH19MWXKFLx48QIZGRkYN24crQdlZGRAXl6+Xn8eY2Nj3LlzR+SxaWlp9XoBSUhIwMDAAGlpabCwsKBlT1NQUlICGRkZoTNPfuJZXV0NSUnJFi0obGNjg4SEBJSWlgqdAU+cOBH379/HpUuXMGbMGL69loyMjBAREQEmk4lOnTrht99+w59//slz7H9nqJx2JpWVlcjOzkZ2djaePHkC4EsuvYODAwDg+fPnCAsLg7u7O3r37t3k1agIBGHQEs8tW7Zw/+3g4IDg4GA8ePAAxsbGGDhwIK0HVVZWNnCKyMrK8lzeCRtbUVHRoLgvv3sJY+/evdi3b5/I1wFf9mUVFRX5zjyFFQTJzs6Gjo5Oi6e4+vr64u+//4a7u7vAwPnu3bujd+/eePDgAc6dO4effvqp3vni4mLExsYiJiYGrq6uCA8Px7hx42BqagoFBQXY2trCxsYGJiYmMDQ0hL6+foNY08ePH4PNZqO8vByfPn1Cfn4+srOz8ebNG7i4uHDHxcTEYN++fdi3bx8UFRXh7u6OIUOGYPDgwd/cA4pAoAPtUCUOxcXFMDQ0hK2trUjXycrKNljiV1ZW8oxvFDZWTk6uQcgUv3sJY968eZg3b169Y9nZ2bSCvXNzc6GsrMxXcDhN3/jteX748AH29vYtXkxYXV0dSkpK+Pz5s9CK81OmTMGDBw8QEhICX19fqKurIzExEZcvX8b9+/e5oVs6OjqQl5fHp0+fMHz4cIwcOZL2e1IUBQUFBSgoKMDIyIg72/waNzc3FBcXIyYmBm/evEF4eDjCw8MhJyeHsWPHYtu2baL/IAgEEaAlniUlJdi5cyfCwsK4szsdHR3MmjWL9rLdyMgIJSUl+PTpE/cDmpaWBiMjI5HHGhsb4/bt29zxLBYLmZmZMDY2pmVLU1FaWso32qC0tBSlpaWQlZWtV43qa2pqapqthqeojBkzBkFBQfDw8BAocmZmZnB1dUVcXBx27NiBgoICZGRkAPhSC9TBwQGenp7o27cvampqcOHCBfz4449N/geiV69e6NWrF9asWYPMzExER0cjJCQEDx8+rOcAq6mpQUlJCdTV1Zv0+QQCLYfRkiVL8PjxY2zbtg1hYWG4dOkSfvrpJ+zZs4f2kldRURHOzs7Ys2cPSktLuXtWQ4YMAQD4+/tzw5GEjfXw8MDDhw8RGxuLsrIyBAYGQl1dHZaWlo35GTQaFovF19HCmXVqa2vzDUNisVh8y9SJGzk5OZiamtLqEDBlyhRQFIXExERkZGRAXV0dU6ZMwdmzZ7Ft2zZ4eXlBXl4eysrKoCgKqampzWq7oaEhpk6disuXLyMhIQELFy7knouIiICDgwMWLFjAdSgSCE0BrZlnQkICTp48iR49enCPWVhYQE1NDatXr6bdPXPTpk1YsWIFXFxcoKysjDlz5nCXZFlZWSgsLKQ1Vk9PD1u3bsXGjRuRm5sLCwsL7N69mytSRUVFGDx4MGpra1FeXs71zO/YsYPeT4UmgrKLhC3ZOYg7OF4QQ4YMwf79+6Gtrd1gH7aurg7Xr1+HhIQEBg0ahMGDByMiIgLm5ub4888/+Tps+vbti7CwMJiYmIjFqWNoaFjv++TkZO4M+MKFCxgyZAjmz59f73eZQGgMtMRTR0eHZyygiYkJrQyVr+9z9OhRnuf+G1gtaCzwpd7k0KFDeZ5TUVER2PGzqaAjnvycRZz40JZ2Fn2NhIQEevfujTdv3tTripmUlIQDBw4gJSUFqqqqcHd3x5QpU3Dr1i28fv0aycnJfMVIUlISRkZGuHnzJgYNGiSuV+ESEBCACRMmIDAwEGfOnEFERAQiIiLQv39/LF26lIgoodHQWrbPnj2bZwfNpKQkWFtbN7lRbQVB4inM056Tk9MqPO3/xcnJCR8/fkR5eTkKCwuxZcsWLFy4ECkpKdDQ0MBPP/0EaWlpaGhoYPz48QCAQ4cOCUxT7d69O54+fYqSkhJxvUY9jIyMsHnzZiQkJGDWrFmQl5fHrVu3EBcX1yL2ENoHfGeeFhYW9fbq2Gx2gz1FNpvNtzVtR0BQ7yJhy/acnBz06tVLaE1NcUNRFHx8fLBt2zbExMSgtLQU0tLSGDduHMaPH18vomHs2LEICwtDSkoKoqOj4eXlxfeeTk5O+Oeff+Dv799i76ylpYXVq1djzpw5OHbsGKZNm8Y9FxcXB3Nzc5HSjQkdG77Kd+LECXHa0Sapqanh60UWtmyvqalptoZv34q2tjY3cN7e3h4LFy7k+R5ycnKYPn06duzYgaCgILi4uPCNsdTU1MTz58+RmpoKMzOz5n4FgaipqWHRokXc70tKSvDTTz+BxWJh0aJFmDp1Kgm6JwiFr3j27t1bnHa0OTj517xmUbW1tcjLywNFUdDS0uJ5PYvFanXB3CwWi7sPGxgYiEOHDmHOnDkCnVoDBw5EaGgoXr58iZMnT2LWrFl8x7q5ueHy5cuYO3duq6o4X1ZWBnt7e0RHR2P9+vU4c+YMNmzYADc3t5Y2jdCKob1+Sk5OxpIlSzBmzBiMGTMGS5YsQXJycnPa1qopKCiAgoICz5lnXl4e6urqoKGhwXdPk81mtxpPe3l5ORYvXoxly5Zxjzk6OmLBggVC894ZDAbmz58PiqJw8eJFgVWapKSkYGdnhwsXLggtIi1OdHR0cOLECRw/fhzGxsZISUnBhAkTMHPmTFqhW4SOCS3xvHbtGsaPHw8Gg4HRo0dj9OjRYDAYmDBhAq5evdrcNrZK8vPz+aZmCqumxGKxwGAwWjSnnUNSUhKGDBmCs2fP4uLFi/WqxdvY2KCiokJoX/Zu3bph+PDhqKurw549ewQ6j/T19VFVVYWnT5822Ts0FQMHDkRMTAyWLVsGOTk5XL16FT/++COp8ETgCS3x3L17N5YtW4bt27dj8uTJmDx5MrZv346AgADs3r27uW1slXDy2nkFwAvb7/zw4QO0tLRaXDzPnTuHESNG4O3btzA3N0d4eHiDOMlRo0bh33//FSiIADB16lQoKyvj2bNnCA8PFzjWxcUFkZGR9eJ6WwsyMjKYN28ebt++DW9vb6xdu5b7/5iIKOFraIlnTk4Ozz1QR0dH7iyro1FZWcnXqSBMPLOysmBoaNhiOe3V1dVYsWIFFi1ahKqqKkycOBHh4eE8K1JpampCU1MTubm5Au+ppKSE+fPnAwACAwMFVpxnMBjo168fTp48ybOtS2tAV1cXgYGB9X7vAwICsHXrVlRUVLSgZYTWAi3xNDU1RXR0dIPjN27cgImJSZMb1RagE+PJb9leXV3dYgWQAWDPnj04fvw4pKWlsXPnTuzYsUNgUZVRo0bh2bNnQoWuX79+cHNzQ0VFBXbu3ClwpqasrAxDQ0NcuXKlTczo3r17h+DgYOzbtw8DBw4USxIGoXVDSzyXLl2KgwcPwt/fH9u3b+cu34OCgsRSrb01wmaz+easC5t5trSnfdasWejXrx8uXryICRMmCB0vJSUFDw8PPH36VKjQzZ8/H0pKSnj8+LHQ/XAzMzMUFxcjISFBJPtbAhMTE1y6dAnm5uZIT0+Hr68v1q1bR2ahHRha4uni4oLr16+jZ8+eyMrKQlZWFuzs7BAREVGvxmJHora2lqeQsNlsoeJZV1cn9lCd+Ph4VFVVAfjSPTM4OBh2dna0r7e1tUVlZaXQLCFVVVVuib9Dhw4JXL5zgufv3bsnsK1xa8HR0RERERFYuHAhGAwGgoKCMGjQIFLhvoNCSzyXLVuGmpoa/PLLL9wCtIsWLeIrDh0Bfu1HioqKUF5ejk6dOvFcmnM87eJKy2Sz2Th06BDGjx+PX3/9tdFLZIqiMH78eNy7d09omFH//v3h6uqK8vJy7Nq1S2iok6enJ4KDgwX2RmotSEtLY8mSJQgNDYWZmRnevXuH06dPt7RZhBaAlnhGR0e3iX0pccIvu+jrtEx+nnhtbW2xeNpra2uxcuVKbNiwAWw2+5sze5SUlGBlZYU3b94IHEdRFObPnw9FRUU8evSIZ5+qr5GWlkb//v1x+PBhWs3oWgM9evRAREQElixZgnXr1nGPi1Ioh9C2oSWezs7OZGnyH/jltQtrN5ydnQ19ff1mT/8rKyvD1KlTcfz4ccjIyODAgQOYN2/eN7c49vT0RGZmplCRU1NT45YqPHDggNAWx8rKyujZsyeOHz9Ou6lgSyMjI4OFCxdyC1pXVlZi6NCh2L59OxHRDgCtqh5du3bF9u3bAYCnYIwaNapprWrlcISD18yTjqe9uXPaP3/+jB9++AGJiYlQVVXF0aNH4ejo2CT3ZjAYGDt2LEJCQtCvXz+BRT48PDzw6NEjREZGYsOGDdi/f79Ar76uri6YTCbOnz/PTcpoS9y5cwevX7/Gq1evEBkZid27d9cr7UdoX9ASz8uXL0NBQYFn1XiKojqceObn50NJSalR2UW1tbXNXj1+165dSExMhL6+Ps6cOdPk4WR6enrQ0tJCVlYWzzYqHDjL95SUFGRkZGD37t0ICAgQOPs1MzPDkydPEB4ejuHDh3/zTFmceHl54eLFi/jll1/w8uVLeHt745dffsGcOXM6dPWx9gqtP+0xMTF8v3jFf7Z3Pn78yDc1U9iync1mN3uY0ooVK/Ddd9/h8uXLzRaHO3LkSLx8+VJoqI6cnBzWrFkDWVlZREVFISIiQui9e/TogcLCQty4caPN7bX37t0bkZGRmDJlCmpqarBjxw6MHDmy2VuREMSPUPHMyMjAsWPHcP78eYFhJx0JzsxTUGomL/Hk9GlvDk97SkoKdztBTk4Ov//+e7NGQ0hKSmLChAm4e/euUO+7kZERFixYAOBLq+d3794JHE9RFBwcHPDhwwfExsY2mc3iQl5eHps3b8aZM2egq6uLJ0+ecHvPE9oPAsXz33//xejRo7F7925s3boVQ4cOJY4j8HcWMZlMFBcXQ0ZGhme3Ro6zqKk97Xfv3oWPjw/mzJkjVmeLjo4OzM3NaTVW8/LywtChQ1FdXY3169ejrKxM4HiKouDo6IjU1FTcuXOnqUwWK+7u7oiOjsamTZvg6+vLPc7pQEto2wgUzz179sDT0xOPHj3Cw4cPMWLECKxevVpctrVa+BXJ+Lr1Br9ZqZaWVpN62m/evIlJkyaByWRCWlpaaAGPpmbgwIHIzc1FcXGx0LFz585Fly5dkJOTg23btgmdsTIYDDg7OyM5OblNZCHxQklJCVOnTuX+PqSkpKB37944ceJEm9uSINRHoHg+e/YM/v7+kJCQgISEBBYtWoS0tDTk5+eLy75WCb+Wwxzx1NPT43ldbW0t3+V+Y7hx4wamTZuGqqoqTJo0Cfv27RN7TySKojB58mTcu3dP6KxXRkYGa9euhaKiIhISEhAUFCT0/gwGA66urnjy5EmbnYF+TWhoKAoLC7F8+XJMnDiR1AttwwgUz4qKCqipqXG/V1JSgpSUlNAlV3uHn0jQqePZVGFK165dw8yZM1FdXY3p06dj27ZtLValSUlJCU5OTnj+/LnQ2ZS+vj7Wrl0LCQkJXLhwAaGhoULvz2Aw4O7ujnfv3rWZQiL8WLx4MQ4cOAAVFRXExsZi4MCBOH/+fJt+p46K0PiJ48eP1/vA19XV4dSpU1BRUeEeo9u3vb1QU1PTqDAlNpvdJDntd+/e5fbc+emnn7B69eoWD+np06cPXr9+jdzcXKGOqp49e2LhwoXYuXMn9u7dCxUVFaEtLxgMBhwcHJCSkoJjx47hhx9+aLE/Ft/KyJEj4ezsjKVLlyIyMhILFy7EtWvXsH37dtKArg0hcObp6OiI5ORk3L9/n/vVq1cvvH79mvv9gwcPxGVrq4HTc/2/CApTYjKZUFBQaBJnUa9evdC3b1/MnTu3VQgn8GX5PnHiRCQmJtJamQwdOhT+/v6oq6vDli1bkJiYKPQaBoMBc3NzGBgYYN++fUIr3LdmOnfujKNHj2LXrl1QVFREbGxsm36fjojAmSevXu0dnaqqKrDZbJ7imZWVBeDL0vS/pKamwtjY+JvEk1MGT1ZWFsePH4eUlFSrEE4OUlJS+OGHH3D69Gl4enoKDQyfNGkSiouLcfnyZaxZswa//fYbz4LMX0NRFPT09KCgoICgoCCMHTu2QfX7tgJFURg3bhz69u2LFy9ecGNy2Ww2SktLW7TmK0E4bSv/rRXAr3dRaWkpioqKICsry3PpVVJSgs6dOzc65fDUqVOYOnUqt6yctLR0qxJODpqamnBzc8OjR4+E7uNRFIWff/4ZHh4eqKioQEBAAF69ekXrOSoqKhgwYABCQ0Pb/OpHT08PgwYN4n5/9uxZ9OvXD1FRUS1oFUEYRDxFJDc3l2dqJmfWqaenx7cdcWOdRceOHUNAQAAiIyMRExPTqHuIE3t7e6iqqtISQgaDgaVLl8LNzQ1lZWVYunQp7a6ssrKycHd3x8uXL1tdR87GwmazcfXqVeTl5cHf3x+LFy8WWkOV0DIQ8RSRT58+8Qw3ys7OBgAYGBjwvK6urq5ROe1BQUFYuXIlAGD9+vUYOnSoyPdoCUaMGIHPnz/T6nElKSmJlStXon///igvL0dAQADtjBxJSUk4Ojpyay+0xqZyokBRFI4dO4bVq1dDRkYGZ8+exYABA3D16lXikW9lEPEUEX6edk4ldF7iWV1dDWlpaZH3Ow8ePMitFbl582bMmDFDdINbCAaDgR9++AFPnz6lFUAvKSmJ5cuXw9PTExUVFVi+fDlu3rxJ61kURcHY2BguLi44duxYm8+Ck5CQwKxZsxAREQE7Ozvk5uZi5syZmDZtWpsoGN1RIOIpIvwyeATNPNPT06Gvrw9ZWVnaz9m9ezc2bdoEiqKwY8cOTJkypVH2tiTS0tKYOXMmbt++TSslUUJCAgEBARg9ejRqamqwefNm/PPPP7Sf16lTJwwcOBDJyck4fPhwm49HNjMzQ0hICDZv3oxOnTohKyurRXtfEepDxFNE+O2rCZp5fvz4Efr6+rTjElksFh49egSKorBz505MnDix8Qa3MIqKipg+fTqioqJotRlmMBj4+eefMXPmTABfWX+5rwAAIABJREFUZt+HDh2inXYqKSkJBwcHdO/eHUeOHEFMTEybXu5KSEhgypQpuHXrFvbv389N7c3Pz0dSUlILW9exIeIpIrxSM1ksFndvj1eYkqjOIgkJCQQFBSE4OBjjx4//NoNbAerq6pgwYQIiIiJoOXU4/ZICAgIgISGB8+fPY8OGDWAymbSeR1EUVFRU0K9fPxQXF2Pfvn1tosGcIDhFWDisW7cOQ4cOxapVq9r8Pm9bhYiniFRXVzeYQebm5qK2thaampo8K6WzWCyhziI2m41z585xy8pxPMntBSMjIwwbNgyRkZG0Z5FeXl7YvHkz5OXlERcXh9mzZwvtn/Q1UlJSsLCwgIuLC27cuIFTp061i1bBdXV10NDQAEVROHr0KNzc3HDixIl2EW3QliDiKQJsNptn47eMjAwA4BmsXVVVBRkZGYHtJ9hsNtavX49FixZh7ty5bXqZKQhra2u4u7vj1q1btAXUwcEBBw8ehKmpKT58+ID58+cjNDRUpJ+RvLw8nJ2dYWRkhP/973+4ePEirS2E1gqDwcD69etx/fp1uLi4cAuNDB06FPfv329p8zoMRDxFoLCwEHJycg287Zzivryqtr99+xbGxsZ8nUVsNhtr1qxBUFAQpKSk4Ofn1yqD35uKXr16oWfPnoiPj6ctgHp6eti7dy+GDx+Ompoa7N69G1u2bBGpLiaDwYCGhgb69+8PNTU1BAUF4cqVK22m2RwvLC0t8ffff+PQoUPQ09PDixcvMHbsWK7zktC8EPEUgffv30NJSanBzJMjnl26dGlwTUFBATQ1NXnW8Kyrq8Py5ctx5MgRSEtLIygoCIMHD24e41sJFEXBxcUF3bp1w82bN2nPQKWlpbFw4UKsWLECsrKyuHnzJmbNmkUrJ/5rJCUloaWlhb59+6JTp04IDAxEeHh4m13yUhSF4cOHIzY2FosWLcK0adO4++5sNptWmBihcRDxFIH3799DRUVFpJknp4bnf6mrq0NAQABOnjwJGRkZHDlyBF5eXs1jeCuDoij0798fzs7OuHr1qkhLaA8PDxw4cAAmJiZ4//49lixZgt9++03kLBxpaWno6OjAxcUFcnJy+OuvvxAWFtZml/NycnJYvHgx1q9fzz0WFRUFJycn7N27t13s9bY2iHiKAKcH0ddUVlYiJycHEhISPPc86+rqeIrn8ePHERwcDFlZWRw7dgwDBgxoNrtbIxRFoVevXvD19UVYWBhtTzrwZW95//79mDp1KqSkpHD9+nVMmzYNN2/eFHm/WEZGBjo6OnB1dYWcnBwCAwMRHBzcLooUcyo1bdu2Da6urjh58iTpJ9+EEPEUAV77Y+np6WCz2TAwMGhQxb2srAxycnI8a3h+//338Pb2xokTJ9qVV11UTExMMHPmTNy8eRO5ubm0r5OSksLEiRMRGBgIW1tbFBUVYfPmzVi5ciXXgScK0tLS0NXVhbu7O0xNTRETE4OjR4+KtLXQ2ti0aRPOnj0LW1tb5ObmYtmyZejbty+OHz/OLTBDaDwS6zj5fwQuJSUlOHHiBPz9/evNGpOSkqCsrFwvzfLBgwdISEhAr169GhT0TU5OhoGBAUxNTcFgMFBRUYHa2lpISUlBUlISPj4+bbacWlOioKCAnj17Ijo6GlVVVdDQ0KB9rbKyMry8vKChoYHnz58jPT0dYWFhyMvLg5mZmcjFpxkMBmRkZKCtrQ01NTWUlZUhJiYGb9++haamZpvL8DEyMsL3338PMzMzpKamIj09HdHR0VBRUYG9vX1Lm9emITNPEaipqWngCRe031lUVITOnTtDQkICpaWlmDRpEmbMmEGWTjxQUFDA9OnTUVVVhbi4OJFmewwGA97e3jh27BhGjBgB4EubEn9/fxw+fLhR+eAMBgMKCgrQ19eHo6MjDA0Nce3aNZw4cQJxcXFtKpyMwWBgxIgRiIqKwqFDh+Ds7Izvv/+ee/7+/fsk0L4REPEUAV4V5AWJJ4vFgrKyMj5//oxx48bh3r17ePXqlUjL046EtLQ0/Pz8YG9vj7CwMHz+/Fmk61VVVTF//nwcPnwYrq6uqKqqwpkzZzBx4kTs3bu3UT93TvFpZWVl2Nvbw9LSEkwmE3/99RfOnDmDBw8etJllPYPBwPDhw3HhwgVu0kZFRQVmzJgBR0dHrFy5EmlpaS1sZduBYrelP6FiIjs7G56enoiOjuaGfRQXF+Pq1auwsrLihh3V1dVh1KhRKC8vx9mzZ+stN1ksFmJjY+Hl5YUff/wRb968gbGxMc6cOUOW6jQoLi7G6dOnoaysDDs7u0bFviYnJyM4OBj37t0D8EU8BgwYgHHjxsHU1PSb7KuurkZ1dTXy8/ORmZkJRUVFdO7cGW5ubiIVgGlpcnJysHTpUty6dQvAlz8WgwYNwtSpU9G3b99GF+/uCBDx5AEv8UxKSkJ2djYMDAy4v1BpaWmYOXMmOnfujODg4Hr3ePXqFQoLC/HXX38hJycHFhYWCA4OhpaWltjfp61SU1OD+Ph4PHz4EH369BFpL/Rr0tLScO7cOcTExHBniVZWVvDx8YG7u/s3t2uuqalBVVUVmEwm3rx5Azk5OSgpKcHV1RXq6urfdG9x8erVKwQFBeHixYvcbSVjY2NcvnyZNKXjAxFPHvASz+vX/197Zx4WVb3G8c8M+w4SiJIyKpdBFhUlQa4iGtcFDUxvZZpZLrndupHZk9zUrKzItBJLr5q3suw+ZpLlkuCOuBBirtcRcnAFGVBAmGFg4Nw/fOY8DkvCiFBwPs9znjlzzu+c8/7mnPme3/q+u3BxccHFxUVMt337dj766COioqJ44403TM7x/fffs2HDBsrKyujbty9fffUVbm5uLZqPtsLNmzdJTk5Gq9Xel9jl5+ezefNmUlJSxNlJzs7O/O1vf2PYsGH3XRqFOzUOnU5HVVUVKpWKmpoaHBwccHV1JSws7A//DGg0GjZs2MDGjRtxdXUlNTVVLPWfOHGCXr16/WmjljY3knjWQ33i+dNPP+Hl5WXi4GPp0qXs2rWL2bNnM3bsWJNzpKamcujQIRwdHfn0009/d267xL0xGAxiT7qjoyNhYWH3DDDXEDqdjr179/LTTz+Rk5Mjbu/evTtRUVFERUU1GD66KQiCQGVlJVVVVVRUVJCdnY0gCDg4OODg4EDfvn3x9va+7+s8CKqrq8nPzxftU6vVDBw4EE9PT2JjY3n88cfp3bt3m55KfC8k8ayH+sQzOTkZb29vE9dyzz33HFevXmXlypX4+/uLjkMqKys5efIksbGxuLu7S2/qZkSv16NSqdi3bx+WlpZERESY3cYoCAIqlYqUlBT27dtnEvrXz8+PiIgIwsLC8PX1bRaRMIppZWUlBoMBtVpNeXk5Dg4O2NnZ4eTkRN++ff+QVf20tDRef/11cnNzxW3du3dnzJgxjBgxgoCAgHYnpJJ41kN94rlp0yYUCoUonjdu3GDixInY29uzZcsWampqSEpK4tq1azz++OMEBQURHh4uNbg/IPR6PZcuXSI1NRWDwcAjjzxyX1XiyspKMjMz2b9/P0eOHDGZzuju7k7//v0JDw8nODi42UICC4KAwWCgsrKSmpoaKisrUavVVFZWYmdnh62tLfb29vj5+aFQKFr9WRIEgV9//ZXk5GS2bt1KYWEhcMdr1alTp8TalTFEdlvHvHpPO0MQBPR6vUk10RgnJyQkhNLSUhYvXszZs2extrYmJyeH6OjoVn/Y2zI2NjaiqBQUFLBr1y5KS0tRKpX1Omi5F9bW1kRERBAREYFeryczM5OjR4+SkZFBUVERO3fuZOfOncCd6aGBgYHi8vDDD5slFjKZDCsrKxOnMW5ublRVVVFVVSWK6/Hjx9m7dy82NjbY2NhgbW2NnZ0dSqWSrl27tthzJpPJCAkJISQkhIULF5Kens727duRy+WicFZUVDB48GBCQkKIjIwkMjKyXgfhbQGp5FkPtUueBQUFHDx4ED8/P/FBX7x4MWlpaYwfP549e/ag0Wjw8PBg3rx56PV6xo0bZ9K5JPFgqa6u5tatW6Snp5Obm4tcLicgIOC+2y4FQSAnJ4djx46RmZmJSqWq4zzE2dmZwMBAAgICCAoKws/Pr8nB/u5lg8FgEBe408N//fp1SkpKsLKyEoXV0tISJycnFAoFnTt3vu+RBE3F+J+4G4VCQWRkJAMHDiQqKupPN0urISTxrIfa4pmWloYgCKL37qqqKsaNG4dWq8XS0hKDwUBgYCCLFi3i7Nmz9OnTh9DQUKnk2QoIgoBWq+XmzZscOXKEvLw85HI5PXr0MKtEWpvKykpycnI4d+4cZ86c4ezZs3Vm58jlcry9venWrRsKhYJu3brRrVs3OnXq1Kzt30ZRra6uxmAwiLOe9Ho9hYWFlJaWIpfLxdKtpaUllpaWuLi44OPjQ6dOnczudPs9Ll68yMGDB0lLSyM9Pd2kLTkjI0PshDpy5AjOzs74+/v/KfsFJPGsh9ri+eOPP9KpUyexp/3QoUPc7RJg9OjRzJ49G0tLS3bv3s24cePabFXlz4QgCJSXl1NcXExWVhaXLl1CJpNhbW1NQEBAs3TMCIJAfn6+KKbnzp1DrVbXO+vIxsYGHx8funTpQufOnenUqZO4dOjQoVlftjU1NaKwGuNuGZsWysvLKSoqMhFXS0tLrKyskMvlyOVyXF1d8fLywsPDAxcXF7PbMA0GAydPnuTgwYOoVCpWrVolnmvo0KGoVCrs7Ozo2bOn2AwSFBSEv7//H36EitTm2Qj0er3JXOaUlBQAfH19efrppxk8eDAAWVlZBAYGmj2YW6J5kclkODo64ujoiLe3NzqdjvLycgoKCjhx4oToSNnS0hIPDw98fX2bXN2WyWSiAD766KPAndLp5cuXUavVqNVqcnNzUavVaDQaLly4wIULF+qcx9raGi8vLzp16kTHjh1xd3fnoYcewt3dXVx3dHRstIjJ5fIGq+yOjo489NBD1NTUiOJqFHuZTCa+dM6cOUNZWRkVFRViqbW2yFpYWGBhYYGTkxNubm64urri7OyMo6MjcrkcS0tL+vXrV8cJSXV1NYGBgWi1Wq5cuUJWVhZZWVni/n/+85+89tprAOTk5JCZmUmPHj3o0aMHHTp0aNRv8KBpUfEsKCggISGBjIwMXFxcmDFjBs8884xZaXfv3k1iYiL5+fn4+/uzZMkS/Pz8gDtvuyVLloixbkaNGsWCBQvq9eZ+L4yzR3Q6HWvXrmXw4MEcPXoUCwsLlixZIpZeqqqq0Gg09zV0RuLBIZPJsLe3x97eHg8PD3r27ElFRQXl5eXcvn0btVpNRkaGOLtGLpdjb29P165dRecujcXa2hpfX198fX1NtpeVlaFWq7l27RrXr18nLy+PvLw88vPzKS4u5vLly78b5dPGxgZ3d3c6dOiAi4sLzs7OJp/GdeN3BweHBkuzRtFr6D/h5OREx44dqampMVlqR481iq1Op+PSpUucP3+eiooK9Hq9eI27RdYounK5nGHDhjF8+HAqKyu5efMm+fn55ObmkpOTQ1BQkHiNgwcPsmDBAvG7q6srXbp04eGHH+bhhx9m0aJF4kulrKwMBweHFuntb9Fq+wsvvICXlxfz5s0jNzeX6dOns2rVKkJCQpqUNj8/n5iYGD766CNCQ0P55ptvSE5OZseOHchkMtatW8fOnTv57LPPAHjppZcYOnQoM2bMaJSdd1fbL168SEpKClu2bKGkpAQXFxdKSkoYOXIkc+fOFY9JTU0lNDSUfv36NWtngUTLYXxRGqdaajQa1Go1hYWFJtVw45/fzc2Nzp074+bmdt9Vbq1WS15eHtevX6ewsJDCwkKKioooKioS15sSs8lop3FAvr29fZ31uz9tbW1NFhsbG/HTOGzK2tr6ni8Ro8gKglBnvbboGjEYDJSVlVFeXo5erxfHSguCgFwuJzs7m5MnT6LRaNBoNCZeyZycnEhMTEQmkyGXy0lISKC4uBhnZ2fc3Nzo0KEDDz30EB4eHkRGRjJgwABsbW0xGAwUFxfj5uZ2z8i2DdFiJc/bt29z6NAhDh8+jJOTE8HBwYwaNYqdO3fWEc97pTWGFzBWl6dOncratWtRqVT4+/uzY8cOZs2aJc4jnzJlCqtWrWq0eBrZunUr69ato6CgAAAPDw80Gg2Ojo5MnjwZuNPmlZ6eTteuXVEqlZJw/okxdqw4Ojri7u6Oj48P/fr1E4cOVVZWirUQnU5HQUEBubm5nDp1ql4XdUaRNQ5JcnJyEpsRHBwcsLGxEUXE3t5erJY2hFarFQW1tLRUXEpKSigpKRHXjZ9arZbbt2+bdNjcL9bW1tjY2Ii/Ve3F2tr6nvuM1X9jabSh0qmFhQUymYyuXbuK41xlMplYWygpKRHHx8rlcgRBoKKiQhx5cevWLdHrGdwZm33t2jWqq6s5ffo0P/zwA56enmzfvt2sURktJp6XLl3C3t4eV1dXcZtCoSAtLa3JadVqtUmHjIWFBV26dEGtVuPv719nv0KhMJkZcS+MwcDef/99BEHAysoKCwsLCgsLsbCw4IUXXkCr1bJ//340Gg2+vr74+Pig0+mkyIXtAKMYODs706NHD5Meb+OnwWAQRdc4372srIy8vDxRfJvq11Umk9Wpjhqr6V26dAEQ2xmNwnO3HXq9XnwB1F6Ms56MM6BqL8Y0xnP9kfjvf/9r8r2hKntqaiqpqakm6TQaDSqVCk9PzyaPPGgx8ayoqKjTFmhra1tvVeReaXU6XZ1Ombv31z7ezs4OnU5X78yHpKQkVq5cWa/Ntb2QG3/cNWvWNJhPCYm2SO3B/G2JmTNnmswmbCwtJp7Gdoa7qaioqHc4wr3S2tnZ1Xn73b2/9vE6nQ5bW9t630gvvvgiL774Yp1z9e7dm5SUlD/l+LP7wdjW2x6R8t5+8+7l5dXk41pMPH18fCgtLaWwsFAsNarVanx8fJqcVqFQcPDgQTF9dXU1ly9fRqFQAHfip2dnZ/OXv/zld6/TEMZSa1OOaUu05zGqUt7bJ+ZMFmixKTBOTk4MGDCAFStWcPv2bU6fPs22bdsYMWIEAJMnT+brr79uVNqhQ4fyyy+/cODAAcrLy1mzZg3u7u707NkTgOHDh7Nu3Try8/O5ceMGn3/+uXishISERHPQouM833nnHRISEoiIiMDFxYU5c+YQGhoKwJUrV0ymuf1eWm9vb9577z3efvttcZznJ598IlbLp0yZwrVr1xg1ahSCIDBy5EimTZvWklmVkJBo40jTMxtAqVSiUqla24wWp73mG6S8S3lvGlLc9t8hLCystU1oFdprvkHKe3vFnLxLJU8JCQkJM5B8pklISEiYgSSeEhISEmYgiWctCgoKmDZtGr169WLQoEHi8Km2TlJSEv7+/gQEBIhL7YigbYmysjLefvttXn31VZPtx48f57HHHiMoKIjRo0eTkZHRShY+OBrK+6RJk+jZs6fJM/D++++3kpXNy+rVqxk5ciR9+vRhyJAhfP755yb7zbrvgoQJ06dPFxYsWCCUlpYKp06dEsLCwoSsrKzWNuuBs2LFCmHu3LmtbUaLsGvXLqFnz56CUqk0ybNWqxXCw8OFzZs3C+Xl5UJycrLQv39/oaysrBWtbV4ayrsgCMIzzzwjbNq0qZUse7CsXr1aOHnypKDX64Xz588LERERwv79+wVBMP++SyXPuzB6c3rllVfqeHOSaDsMGzaMc+fOMWfOHJPtR48excPDg3HjxmFvb8+YMWPw9PTkyJEjrWRp89NQ3ts6M2bMoFevXlhbW6NUKgkNDeX8+fOA+fddEs+7aMibU1M8Mv2Z2bFjB0FBQURGRvLOO+802evPn53a3rigfd1/gDfffJPevXszYsQIvv3229Y254FQXV3NmTNnRPd/5t53KQzHXTTF81NbY8KECTz77LM4OjqSk5PD/PnzWbZsGfPnz29t01oMnU5Xx1FNe7n/AImJibi5uWFhYcGxY8d45ZVX8PT0FMOLtBU+/PBDOnToQFRUFGD+fZdKnnfRFM9PbQ13d3dcXFywsLBAqVQybdo0E+cr7YF7eetq63Tu3Bk7Ozusra0ZNGgQMTExbe4ZWLNmDbt37+bTTz8VnYGYe98l8byLu705GWmqR6a2glarbXeB7BQKBTk5OSbb1Gq16K2rvVGf39w/M8uXL2fr1q18/fXXeHp6itvNve+SeN7Fvbw5tWU++OADTp06hV6v5/z586xdu5bY2NjWNqtFCQsLo6ioiM2bN6PVaklOTqagoIABAwa0tmkPnIKCApYvX45araayspK9e/eye/duYmJiWtu0ZmHu3LmcOHGCjRs3iuF5jJh736XpmbXIy8sjISGBzMxMXFxcmD59uhivqC3z7rvvkpqaSlFREV5eXjz99NM899xzLRKFsKXZu3cv8+fPp6KigpqaGuzt7YmPj2f8+PFkZGTw1ltvkZubS9euXVm4cCHh4eGtbXKz0VDeY2JimDdvHqdPn0ar1dK9e3fi4+MZNGhQa5vcLCiVynodm587dw7ArPsuiaeEhISEGUjVdgkJCQkzkMRTQkJCwgwk8ZSQkJAwA0k8JSQkJMxAEk8JCQkJM5DEU0JCQsIMJPGUaFUmTZpEUlJSa5vR7Fy9ehWlUsnVq1d/N53BYGDs2LF1fGuac677RRAEZs6cyaxZsx7oddoKkni2U8rKyvjwww+Jjo4mODiYwYMH8/zzz7Nly5b7Ou+WLVtQKpXNZOUdJ81Dhw5tUvpJkyY12/UfNF988QVlZWUsWbKktU1BJpOxdOlSfv31V3bs2NHa5vzhkbwqtVNmzZpFYWEhCQkJdOvWjaKiItLS0ti0adN9eZAfPnw4/fv3bzY7J0+ezBNPPNFs5/sjodfrWb9+PQsXLsTGxqa1zQHuTFH+xz/+werVq9vM1MwHhVTybIcUFRWRkZHB3LlzGTp0KN26dSM0NJT4+HjWr18vptPpdCxcuJCwsDDCwsJ45ZVXuHXrlrhfqVTyxRdfMHXqVEJCQvjkk0/YtWsXzz77rJjmm2++4dFHHyU4OJjw8HDi4+O5efNmo2398ssvmTdvnvh90qRJvPHGG7z22msMGDCAIUOGkJKSAsCxY8dYuXIlGRkZKJVKlEolx44dAyA9PZ0xY8bQq1cvYmJiTBxcJyUlMX78eD777DOio6MJCgriyy+/rFPivX37NsHBwfz6668ATJkyhQEDBhAcHEx0dDRffvllo/MFsGfPHmQyGcOGDTPZbjAYWL58OY888gj9+/endnTwrKwsRo8eTUhICH379mXixImcOXMGgAsXLqBUKrl06ZLJMWPHjmXt2rUAHDlyhCeeeEIMSTF//nx0Op2YNi4ujtzcXE6dOtWk/LQ3JPFshzg7O2Nra8vhw4eprq422Wdvby+uL1q0iKtXr7JmzRrWrVuHRqNh8eLFJum/+OILnnzySZKTk+stITo5OTF37lySk5NZuXIl2dnZJCYm3pf927dvJyAggA0bNjBixAgSEhLQarX06dOHZ599lt69e7Nnzx727NlDnz59+O2333jxxReZMGECW7Zs4fnnn2fevHlcuHBBPOeJEyfIy8tj5cqVbN68mdGjR3Pjxg2ysrLENCkpKXTq1Ik+ffoAEBwczMcff8zWrVuZMWMGiYmJHD16tNH5OHr0KNHR0cjlpn/Djz/+mO3bt7Ns2TK++eabOnOsLSwsePrpp9m4cSMbN27E3d2dl156iZqaGvz8/AgMDOSHH34Q02dnZ3P+/Hni4uIoLi5m9uzZxMTE8P3337NgwQKKi4spLy8X0zs6OhIREdGkvLRHpGp7O8TKyoqEhATeeustduzYQa9evQgICGDgwIGEhoYCdzopfv75Z9LT03FycgLg1VdfZcKECVRXV4tOFt57773f9T5T2zPTU089xYYNG+7L/ueff57nnnsOgOnTp7N+/Xpyc3MJCAjA2dkZGxsbE8/g69atY8KECTz55JMA+Pr6smvXLlJSUvDz8wOgX79+vP322ybX+etf/8q2bdvo27cvcEe0785PfHy8uN69e3e+/fZbTpw40WhHIjk5OYwePdpkm16v56uvviIpKYnIyEjgjr/JpUuXiml69+5N7969xe9z5swhNjYWjUZDx44dGTt2LOvXr+ell15CJpOxZcsWBg0ahKenJ6dOnUKr1TJy5Ei8vLzo0aNHvW3Kvr6+ddy0SZgiiWc75amnniIqKop9+/Zx7tw5Dhw4wKpVq4iNjWXp0qVkZ2ej1+uJiIgQjxEEAYPBgEajwcvLC6BOqak2V69eZcOGDWRkZKDRaCgpKcHDw+O+bL/b05Obmxtwp0rdECqVCpVKZVKtNhgMJq7J6vO4Exsby5IlS0hISKCkpIRjx46ZlLz37NnD9u3bOX36NKWlpZSWljJkyJBG56O4uLiOe7TLly+j1+tFwa6P8vJyvv76aw4dOsSlS5fEvBsd+o4ePZrExESOHTvGI488wk8//cTChQsB8Pf3x9/fn7i4OAYNGkRISAgjR46kQ4cOJtfo2LEjKpWq0Xlpj0ji2Y7p2LEj48ePF78nJyfz+uuvM3PmTKqrq7G3t2fz5s11jnN3d2/U+XU6HePHj8fHx4cXXngBb29v0tPT+e6775otD41xmVdTU8OUKVMYM2aMyXZjibohoqOjWbhwIYcPH+bKlSv06tWLLl26AJCamkp8fDzPPPMMixYtws3NjUWLFjXJdkEQ6oi2Xq8HEL2c18err76KSqVixowZ+Pr6otVqmTZtmrjf1dWVIUOGkJycjF6vx2AwiKJubW3Nd999x+7duzl69Cj/+c9/+PTTT/nxxx9NHB8LgtAm3RE2J5J4tkNKSkqoqakRS21GfH19gTslGOOfsrKykp49e5p1nezsbDQaDT///DOOjo4AD7wqaGNjUyekQo8ePfjtt9/EgF+NxdbWlmHDhrFt2zauXLlgnk31AAADOklEQVRCXFycuC89PZ2oqChef/11cVtTw3W4ublx48YNk23e3t7IZDJycnIIDg6u97jDhw+zfPlyMbZQfeM/x44dS3x8PLdu3SI2NhYrKyvgTonb2tqamJgYYmJi0Ov1hIeHk5mZaeL0+8aNGyaBECXqInUYtUOuX79OTEwMSUlJZGVlcfHiRfbt28cbb7zBwIED8ff3R6FQMHz4cF5++WX27t3LxYsX2b17N7Nnz270dby9vbGysmLz5s389ttvbNu2jTVr1jzAnN1pezx79iwHDhzgxIkTaDQapk+fzv79+/nggw/43//+x9mzZ8WOoXsRFxdHSkoKZ8+eZeTIkeJ2hULB6dOnyczM5MyZM6xYsYKTJ082ydb62hXd3NwYPHgwiYmJnD59mvT09DpB+BQKBT///DM5OTkcOnSoTm88wKBBg3BwcODAgQMmQ8+OHz/OtGnTOHToELm5uezcuZOqqiqx7dfIhQsXxJepRP1I4tkO6dq1K1OmTCEtLY05c+bw2GOPsXjxYqKiolixYoWY7t133yUiIoIFCxYQFxfHhx9+iL+/f6Ov4+7uzpIlS1i/fj1///vf+eGHH5rUJmgOQ4cOZdSoUbz88svMnDmTGzdu4O/vz+rVq/nll1948sknmTp1KidPniQgIOCe5wsLC8PZ2ZnIyEhcXFzE7RMmTKB///5Mnz6dWbNmodfrmzw5IDw8nNTUVGr7I3/zzTeRyWRMnDiRpUuX1vnN3n33XS5cuMDYsWNZtmyZ2LF0NxYWFkRHRxMYGGhyz7p3746rqyv/+te/iIuLY/369Sxfvpzu3buLaUpLSzly5Eib8qD/IJA8yUtItBJ6vZ4hQ4bw3nvvMXjw4GY//2OPPcb48eOZOHFik477/PPP2bp1Kz/++GOz29SWkEqeEhKthI2NDVOnTmXZsmV1xtveL8ePHyc3N7fOUKh7cevWLf79739L89sbgSSeEhKtyOTJk7G0tGz2ue2bNm0iOjrapKnhXlRXV/Pyyy/Tr18/k/ZdifqRqu0SEhISZiCVPCUkJCTMQBJPCQkJCTOQxFNCQkLCDCTxlJCQkDADSTwlJCQkzEASTwkJCQkz+D+5GygI7aRoagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 324x252 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(4.5,3.5); \n",
    "\n",
    "clrs = ['k','C7']\n",
    "\n",
    "disease = 'COVID-19'\n",
    "Ymedian = [np.median(yy[j]) for j in range(len(x)-1)]\n",
    "Ylower = [np.percentile(yy[j],[2.5])[0] for j in range(len(x))]\n",
    "Yupper = [np.percentile(yy[j],[97.5])[0] for j in range(len(x))]\n",
    "ax.plot(x[1:], Ymedian, lw=2, color=clrs[0], label=disease)\n",
    "ax.plot(x, Yupper, lw=.3, color=clrs[0])\n",
    "ax.plot(x, Ylower, lw=.3, color=clrs[0])\n",
    "ax.fill_between(x, Ylower, Yupper, color=clrs[1], alpha=.35)\n",
    "\n",
    "ax.set_xlabel('Serial interval (days)'); ax.set_ylabel('Probability density')\n",
    "\n",
    "disease = 'SARS'\n",
    "res = fsolve(lambda x: find_params_Weibull(x, mean_SARS, sd_SARS), (1,1))\n",
    "res_ = np.exp(res)\n",
    "Y_SARS = [weibull_cdf(xx+xstep, res_[0], res_[1])-weibull_cdf(xx, res_[0], res_[1]) for xx in x]\n",
    "ax.plot(x, Y_SARS, label=disease, lw=2, color=clrs[0], ls='dashed')\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "ax.set_ylim(bottom=0, top=.003)\n",
    "ax.set_xlim(0, 20)\n",
    "ytks = [0,.001,.002,.003]; ax.set_yticks(ytks); \n",
    "ax.spines['left'].set_bounds(0,.003)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False);\n",
    "\n",
    "my_dpi = 300\n",
    "plt.savefig(\"../../results/Fig-MAIN.pdf\",\n",
    "            format='pdf',\n",
    "            figsize=(3.54331/my_dpi, 3.5/my_dpi), dpi=my_dpi,\n",
    "            bbox_inches='tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Location of the max and mean of the distribution\n",
    "\n",
    "(for explanatory purposes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[31mMaximum (x): 2.93\u001b[0m\n",
      "\u001b[32mMean (--): 4.59\u001b[0m\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAD4CAYAAABsdWSLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd1hU1/b3v2eGNgy9SO8KSFMERRFQQWLvJd5rjMaSaGxJfteWeGOKxpSbYov3YkxMNJqrUaOiWAAbiGJHQVF0KEORIm1gmMZ5/+CduSJTzox09ud55gHmrHPOOuh82XuvtdeiaJqmQSAQCAStYHW0AwQCgdAVIeJJIBAIOkDEk0AgEHSAiCeBQCDoABFPAoFA0AEingQCgaAD7SqepaWlWLhwIYKCghAZGYl9+/bpbJuYmIjY2FgEBgZixowZePTokeLYmTNnMHXqVISEhCA0NBQLFixAXl5emz0XgUDoebSreK5fvx6Ojo5ITU3Fjz/+iO3bt+P27dta25aUlGD16tVYv349rl69itjYWKxcuRLylFUul4t//vOfSElJwaVLl2BnZ4d//etf7facBAKh+9Nu4llbW4uUlBR88MEHMDU1RWBgIMaNG4eEhAStbRMTExEWFoZhw4aBy+ViwYIFKC8vR3Z2NgAgIiICwcHB4HA4qKurQ1VVFQIDA9vrUQkEQg+g3cQzLy8PxsbGsLCwULzn7u6O3NxcrW15PB6cnZ0Vx9hsNlxcXMDj8ZpdZ86cOYiIiIC5uTkWLlzI2FepVAo+nw+pVMr4HAKB0LNoN/FsaGiAkZFRs/eMjIxQX1+vta1QKASHw9F4rb179yItLQ1CoRCffvqpUr+2bdsGHx+fZi9/f3/ExMSgpKRE6+ckEAg9g3YTTyMjoxYjuYaGhhYiyMSWw+FAIpEwupaVlRUWLFiA06dPK/Vr+fLlyM7ObvZKSkrS6tm6IjcKM3CjMKOj3SAQuiztJp5ubm6oqalBeXm54j0ejwc3Nzetbd3d3ZGTk6M4JpPJkJ+fD3d3d6X3rqurg4GBQSs9SfcgPjsR8dmJHe0GgdBlaTfxNDU1xZAhQ7B161bU1tbi3r17iI+Px+jRowEAc+fOVaQjabKNjo7G9evXcfHiRdTV1SEuLg7W1tbo27cvAOCzzz7DlStXUFdXh/z8fPzwww8YP358ez0qgUDoAei15802btyIDz/8EOHh4TA3N8fSpUsRGhoKACgoKEBlZSUjWycnJ2zevBmff/45SkpK4Ovriy1btoCiKACAtbU1Pv30U5SUlMDS0hLjx4/H8uXL2/NRCQRCN4ci9TxbwufzERMTg6SkpGZR/e7EJ8nfNX2N/qCDPSEQuiZkeyaBQCDoQLtO2wmdh2WD53W0CwRCl4aIZw/Fxtiqo10gELo0ZNreQ7mSfwNX8m90tBsEQpeFjDx7KGdzLgEAwl1DO9gTAqFrQkaeBAKBoANEPAkEAkEHiHgSCASCDhDxJBC6MGVlZdi0aRNGjhyJgIAAhIeHY/bs2di3bx/EYjEA4OrVq1i4cCEGDhyIwMBATJgwAXFxcYriOgKBAP369cNvv/2m9B5fffWVYntzdHR0swpkfD4fALB27VrFewEBAYiIiMCiRYsQHx+v1fPQNI0nT55g1KhR2L17d4vjQqEQ3333HWJiYhAcHIxZs2YhLS1Nq3u0FiRg1EP5YOjbHe0C4RXh8/mYNWsWnJycsGrVKnh5eUEgECAlJQXbt29HSEgIHj9+jLVr12L27NlYsWIFuFwubt26hR9++AHp6emIi4uDiYkJRo4ciWPHjuHNN99sdo/GxkacPHlS8f7+/fuxe/duXLp0Cbt374a9vb3CdvDgwdi0aROkUimKi4tx7tw5rFu3Dmlpadi0aROjZ/Lz80NjY6PK48uXL0dFRQU2btwIa2trnD17FgsXLsSvv/6q2L7dbtCEFhQUFNDe3t50QUFBR7tCIKjknXfeoSdPnkyLRKIWxyoqKujCwkI6NDSU3rRpU4vjGRkZtI+PD3306FGapmn68uXLtLe3N52Tk9PMLjU1le7bty9dUlKieG/r1q30uHHjmtmtWbOGfvvtt1vc5/bt27S/v7/iPprIycmhc3Jy6IiICPqnn35qduzp06e0t7c3fffu3Wbvr1ixgl68eDGj67cmZNreQ7nAS8MFXsdMdwivTnV1NS5evIj58+crLbdoZWWF+/fvo6amRmkXhcDAQAwZMkRR5zY8PBx2dnY4fvx4M7tjx44pjulC//79MX78ePzxxx+M7L28vODl5QV9ff0Wx6qrqwE0VV17kbCwMGRlZenk36tApu09FLlwDvcY0sGedD7mzJmD5OTkdr1ndHQ09u7dy9g+Ly8PjY2N8PPzU2mTm5sLExMT9OrVS+lxLy8vXLlyBQDAYrEwceJEHD9+HO+99x4oioJQKMTZs2exceNG7R7mJQICAnDu3LlXugYAeHt7w8zMDDt27MC6detgZWWFZ8+eISsrS2lHiraGjDwJhC4IzbAYmrxMI5PjU6ZMQVFREdLT0wE0NVpksViIiYnR3dGX7lNYWIjAwMBmr8LCQkbXMDY2xo8//ogHDx4gPDwcfn5+mDp1Kh49etRiNNoekJFnJ0YgEODmzZt4/vw5JBIJbG1tERUVBTab3dGudWu0GQF2FC4uLqAoCvfv34eXl5dSGzc3N9TW1qK0tFTp6PPJkyfNOjl4eXkhKCgIx44dQ1hYGI4dO4YxY8a06CemLVlZWfDx8QEA9OrVC3/99Vez46pGxsoYOHAgTp48ibKyMshkMvTq1QurVq3SeVnhVSAjz05IaWkpDh48iIMHDwIAHB0d4erqitraWnz99dfIz8/vYA8JHY2VlRXCwsKwY8cONDQ0tDheU1OD8PBwmJqaIi4ursXxjIwMpKWlKbozyJk8eTLOnDkDPp+PK1euYPLkya/k5/3793HixAnMnj0bAKCvr69Y11S3vqkJW1tb2NvbIzc3F2fPnsWkSZNeyU9dICPPTsb169eRl5cHZ2dnUBQFU1NTxUjT3NwcTk5OOHToECZMmABvb+8O9pbQkXz88cf4+9//jmnTpmHp0qXw9vaGUCjEtWvX8Ouvv+LAgQNYv3491q1bB6FQiEmTJsHU1FSRqhQVFdWiPc24ceOwefNm/OMf/4CTk1OL9J+SkhLU1NRAIpGAz+fD3t4eenpNMtLQ0AA+n4+GhgaUlZUhKSkJf/75J6ZPn46xY8cyeqaamhoATSlSIpEINTU1MDQ0hKGhIQCguLgYZWVlMDAwQEZGBrZs2YLRo0dj5MiRr/rr1BpSSV4JHVVJvrq6GqdPn4aLiwvMzMxUTs8lEglOnjyJJUuWNOttrw0iaVMCtaEeaYzXlSksLMTOnTuRkpKCsrIycDgcBAQEYMKECZg0aRL09PSQlpaGuLg4ZGRkQCQSwdXVFZMnT8Zbb72ldNS3YsUKnDlzBsuXL8eyZcuaHYuOjm62Rin/jKxduxZHjx4FAOjp6cHc3BxBQUGYOXMmoqOjGT+PfHr/IsuWLVO00UlLS8OiRYugr68PT09PzJgxAzNnzgSL1f6TaCKeSugI8WxsbMS+ffvg4+MDLpercV2zuroaKSkpWLZs2SuvSREIBO0ha56dhPj4eLi6ukJfX59RQMjc3Bx9+/bFsWPH1O7IUMWZxxdx5vFFXVwlEAgg4tkpyMrKApvNBpfLBYfDYXyep6cnioqKUFpaqvU90wpuIq3gptbnEQiEJoh4dgLu378Pa2trcLlcrc8dMmQI4uPjGef9EQiE1oGIZwdTUFAAPT09GBoaakxoVoa5uTmEQiEKCgrawDsCgaAKIp4dTGpqKpycnBSpGLowdOhQHD9+XFFijEAgtD1EPDuQzMxMGBkZQU9PT6dRpxwOhwOZTKbT2ieBQNANkiTfgTx8+BCOjo5aBYlUERYWhoSEBLz11luMovWfRH/wyvckEHoyZOTZQVRVVUEmk8HAwOCVRp1yzM3NIRAIFGW7CARC20LEs4NISUmBq6trq4w65fTp0wcpKSmMIu/HH57D8YevXiaMQOipEPHsIGpra8FisVpl1CnH09MTT58+RV1dnUbbW0X3cKvoXqvdm0DoaTAST3kjKULr8PjxY5ibmyutAP6qmJiY4NGjR61+XQKB0BxG4hkREYGNGzd2SKn77siDBw9gY2PzSulJqggLC0NKSorSMmUEAqH1YCSeX331FcrKyjBr1ixMmjQJv/32GyorK9vat26LQCBos2vr6emBpmny70MgtDGMxHPEiBHYsmULUlJSMHPmTJw4cQJRUVFYsWIFLl68qFNhip5KXl4eDA0NWzVQ9DK+vr5IS0tTGzgyYOvDgK19EVoCgdCETiXpCgoK8N133yEhIQFAU1XnyZMnY9q0aXB3d29tH9udtixJd/z4cdjZ2cHMzKxVr/sy586dw/z582FiYtKm9yEQeiqMo+319fU4cuQI5syZg1GjRoHP5+Pjjz/GlStX8NFHHyEzM5NxteiejEgkatUIuyooikJFRUWb34dA6Kkw2mG0evVqnDt3DlwuFxMnTsSGDRvQu3dvxfHRo0dj9OjRpDiFBiQSCUQikaJtQVvi7++PtLQ0ODs7K91x9GfmKQDAdH/yB49A0AVGn+KGhgZ8//33iIqKUlvu3sXFpdUc645kZma2WZT9ZRwcHHD//n3U1tYqbdVx/9lDAEQ8CQRdYTRt9/b2RmhoaAvhzMnJwYEDB9rEse5Ifn4+rKys2q3fCpvNRllZWbvci0DoaTD6FO/YsUPR1e5FJBIJvvnmm1Z3qrvS3rmXAQEBSE9Ph1Qqbdf7Egg9AUbiSdN0iyCHVCpFSkoKjI2N28Sx7ohAINCpR7Wu9OrVCxUVFW2aV0og9FTUrnn6+vqCoihQFKW0fSibzca6devazLnuREVFBfT09NpkS6Y62Gw2KioqWqx7mhhq3/KDQCD8D7Xi+dtvv4GmacydOxffffcdbGxsFMcMDQ3h6uoKS0vLNneyO3D79m14eHi0e3/poKAgXLlyBe7u7s2i7v8Y+k67+kEgdDfUiuegQYMANBXtJbwatbW1bZ4YrwwbGxvcunULAoEA5ubm7X5/AqG7olY8t2/fjnnz5iExMVHtRSZPntyqTnVHGhoa2iU5Xhl6enooKytrJp77M/4CAPw9iPzbEQi6oFY8jxw5ghkzZmDr1q0qbSiKIuKpgcbGRtTV1bX7eqec4OBgpKamwtPTU7Fs8Kj8aYf4QiB0F9SKZ3JycrOvBN14/PgxevXq1a6R9hextLREVVUVBAJBhywdEAjdEUbRi4qKCmRmZip+zsnJwS+//IJr1661mWPdiUePHsHe3r7Dpu1A09Sd7HUnEFoPRuL5+eefY9++fQCAkpISzJgxA3/88QeWLl2KQ4cOtamD3YHO0E89KChIY5k6AoHAHEbieevWLYwfPx4AcPr0abi4uODMmTP4+uuv8fPPP7epg90BoVDY7ilKL2NjY4OKigoIhUIAgJWxJayMSZoZgaArjAqD1NXVwdHREUBTvmJYWBiApoZjxcXFbeddN0AgEICmad3XO2kaeHG6//LPWl2KRm1tLYyNjbFi8Fu6+UMgEAAwHHn6+/vjzz//xK1bt3D58mVEREQAAJ48eQJbW9s2dbCrk5GRAWdnZ53E02b7dth9+WWTYAIATcPuyy9hs327Tr54enri5s2bOp1LIBCaw0g8161bh/j4ePz9739HaGgooqKiAABxcXEKISUop6KiAlwuV/tgEU2DXVsLq717Ybt5M+rr6tBr82ZY7d0Ldm3t/wRVC7y8vPDkyRM0NDRgz62D2HProNbXIBAITTCatvft2xcXLlxAZWUlrKysADTlLn733XeKnwnKkUqlugVpKAp5772HzMxMRO7bB5v/H7BLCwuD0XvvwegVIvc1NTXIreLrfD6BQNCiDQdFUc2EksViwcnJqU0bmXUHdA0WFRcXY9ny5Zj10pry9Px8vLt0qc5V+21tbUlfdwKhFWD0qX7y5Anmz5+PgQMHom/fvi1eTCktLcXChQsRFBSEyMhIRfqTLraJiYmIjY1FYGAgZsyY0UwQkpOTMX36dISEhGDw4MH48MMPUV9fz9jP1kIoFEIqlWrddqOmpgbr1q1DLo+Hf73UmfQ7mkZ+Xh5WrlwJHo+ntU8BAQG4c+cOSVkiEF4RRp/qVatWwdzcHJ9//jksLS11TvZev349HB0dkZqaitzcXCxatAj+/v4IDg7WyrakpASrV6/G999/j9DQUPz+++9YuXIlTp06BYqiUFdXhyVLliAsLAwNDQ14//33sWPHDqxatUonv3UlOzsbdnZ2WgWLaJrGt99+C35BAb6nKMx89gzP58zBs7VrYffll3h9714Ye3picXU1Pv74Y/z4448wNTVlfH02mw2apiGVyaDfDr2UCITuCqNPT35+Pvbv3w9vb2+db1RbW4uUlBRcuXIFpqamCAwMxLhx45CQkNBCPDXZJiYmIiwsDMOGDQMALFiwALt27UJ2djZ8fX0xYcIExbVMTEwwevToDtliWlhYCBsbG62m7VeuXEFqaiqMuVwMGjoUzxsb8WztWoCimr4CGGhsDO/sbDx69AjffvstNmzYoNUfNC6XC24jDWtTa62fiUAgNMFIPAMDA5Gfn/9K4pmXlwdjY+NmRXnd3d1x+fJlrW15PF6zfupsNhsuLi7g8Xjw9fVtcb27d+/C09NTqV/btm3Ddh1TfzQhlUq1EjWZTIZdu3YBAN566y0Ip0yB8MW8TrmAUhTWFxVh8eLFSElJQXJyMmJiYhjfJzg4GNeuXcPC2FlaPQ+BQPgfjMRz6NCh+Prrr1XWgxw4cKDGazQ0NMDIyKjZe0ZGRkrXIjXZCoXCZoWZ1V0rOTkZFy9exNGjR5X6tXz5cixfvrzZe3w+XysxUoW22zKvXLkCPp8Pe3v7/42eXxbf//+zo6MjlixZgm+//RY7d+5EaGgo43qdRkZGEIvFqKurI4VCCAQdYSSe//rXvwAAc+bMaXGMoig8ePBA4zWMjIxaNCJraGhQGq3XZMvhcFoIk7JrXb16FevWrcO2bdsUO6Tak4aGBqU901Vx/PhxAMC0adMYBZlGjx6NxMRE3L17F7t378YHH3zA+F5PLZ7h3+n7sHrku4zPIRAI/4OReLZGJXk3NzfU1NSgvLxcMWrk8Xhwc3PT2tbd3R2XLl1S2MtkMuTn58Pd3V3xXmJiIj7++GPs2LEDoaGhr+y/tgiFQshkMsbBovz8fNy+fRtGRkaIjY1ldA5FUVi5ciUWLlyI06dPY+rUqc1+B+pgmRvgSSlPaXM/AoGgGa0SEGmaxrNnzyCTybS+kampKYYMGYKtW7eitrYW9+7dQ3x8PEaPHg0AmDt3riIdSZNtdHQ0rl+/josXL6Kurg5xcXGwtrZWpE398ccf+OKLL/Dzzz93iHACTWJvY2PDOE3pzJkzAIARI0bAxMSE8X1cXV0xfvx4NDY2Ii4ujvF5emw2ZLKmIs0EAkF7GIlnVVUV3n//fQQGBmLEiBHIz88HACxdulSrYMvGjRtRUFCA8PBwLFmyBEuXLlWIW0FBASorKxnZOjk5YfPmzfj8888RFhaGpKQkbNmyRTGCOnnyJIqLizF16lT4+fkpXunp6Yx9fVXy8vJgZWXFaFRH0zQuXLgAAIxHnS8yZ84cGBsbIz09Hffu3WN8HkU15dMSCATtYTQs+vzzz1FRUYHff/8db7zxhuL9KVOm4LvvvsOyZcsY3czBwQG//PKL0mMvpxKpswWAMWPGYMyYMUqP7d27l5E/bYlYLGa83pmVlYVnz57B1tYWAQEBWt/L0tISU6dOxb59+7B//35s3ryZ0XkmXBOkpqbCw8ODTN0JBC1hNPJMSUnBmjVr0K9fv2YfMi8vLxQWFraZc12ZlwNe6rh48SIAYPjw4TrX/ZwyZQqMjIxw/fp1ZGdna7S359jCxcwRlZWVHbL7ikDo6jD+pDa+tE0QAHJzc7Xa3dKTEIlEjEdz8nYmQ4cO1fl+5ubmmDhxIgBg//79Gu3HOEZhjGMUKIpqtlxCIBCYwUg8R40ahR9++AECgQBAU5T34cOH+OqrrzBq1Kg2dbArQtM0RCIRo0g7n89HYWEhTE1NtaoToIzp06dDX18fqampePqUWXdMf39/XLly5ZXuSyD0RBiJ54cffghra2sMHToUYrEY06ZNw5QpU+Dj44P/+7//a2sfuxylpaUwNjZmFGm/fv06ACA0NFSrnFBlWFlZYezYsQCaMg7UcTj/DA7nn4GDgwNKSkrQ0NDwSvcmEHoajAJGRkZG+Prrr7FixQrk5ORAJpOhd+/eSnM0CU3dRXv16sVo/VI+ZR80aFCr3HvmzJk4ceIELl68iLfffrvFTiw5NRJB859ralrs6iIQCKrRKjrh7OyM4cOHIyYmhginGqqqqhhVjxcKhbh79y4oimq1fFQ7OztERERAJpMpdixponfv3u2axkUgdAdUjjznzJnDOODx22+/tZpD3QGmkfYHDx5AIpHA29sblpat18ly6tSpuHTpEuLj4zF79mwYGhqqtXd3d0diYiJEIpFGWwKB0ITKkWdYWBgGDRqEQYMGQSAQQF9fX/Gz/PX48WM4OTm1p79dAmWZCcqQb3v18/Nr1fv7+/vD29sbNTU1SEpK0mgv/yNZU1PTqn4QCN0ZlSPPFxPfT5w4gQ0bNqBfv37NbGxsbEikVglMqynJC6q8apT9ZSiKwrRp07B582YcOXIEY8aMaTGLcDa2b/azq6sr7ty5o9MOJwKhJ8JozbO4uFjpdG7AgAGKgAfhf4hEIo2RdpqmFeKprAbpqxIVFQVra2vk5ubi1q1bLY7HOgxFrMP/8kr79OmDhw8fQiwWt7ovBEJ3hJF4BgcHY8+ePS0KgiQnJ8PY2LhNHOuqiMViSKVSjWlHz549Q1VVFczMzNqkXJ6+vr4iaf7IkSMa7VksFmiaRnV1dav7QiB0RxilKm3cuBFz585FTEwMAgMDoaenh5ycHOTl5eHrr79uax+7FEVFRbCwsNA48nxx1NlW+8rHjRuH33//HdeuXQOfz29Wff+P3JMAgFnu4xTvOTg4IDMzE8OHD28TfwiE7gSjkaeLiwtOnz6NlStXwsXFBdbW1pg2bRoSExMVZeIITeTl5cHa2lpjjqc8WNQWU3Y5FhYWior4L1fSF8oaIJQ1T4z38/NDRkYGRCJRm/lEIHQXGLdPNDAwwJQpU9rSl26BQCCAra2tRrv2EE+gqWBIQkICzpw5g7feekttrVD51L2mpobRMxAIPRndSvgQVMKk6ZtEIlH0mW9r8fT09ERwcDAaGhoUBZfV4ejoqFVNUAKhp0LEs5VhUmX/6dOnkEgkcHZ2bpcGbJMnTwbQ1CNJUw6qn58f7t27R6buBIIGiHi2MlKpFDRNq7Vprym7nMGDB8PW1haFhYW4ffs2AMDDxAUeJi4tbFksFiiKIlF3AkEDRDxbmYaGBsaR9tZOjlcFm83GuHFNUXX5fvfhdoMw3E55MRInJycydScQNMBIPHfv3o2Kioq29qXLI6/jqUk85eudPj4+7eEWAGDs2LHQ09NDWlqaxr5Fvr6+uH//PkmYJxDUwEg8T506hWHDhmHx4sVITEzUqXtmT6C6uhpGRkZq05REIhH4fD5YLBY8PT3bzTcrKytERkaisbERJ0+exF7eMezlHVNqK/e/qqqq3fwjELoajMTz8OHDOHr0KDw9PfHpp58iMjISX375pWIERWiisLAQlpaWancX5efno7GxEc7OzjAwMGhH76DYcXTq1CmIZRJIG1VXf+rTpw/S0tLayzUCocvBeM2zT58+WL16NS5evIjNmzeDx+Nh0qRJmD59Og4cOIDa2tq29LNLUFRUBCsrK7U28vYYHh4e7eFSMwICAuDu7o7KykqNASEPDw/k5uaS5nAEggq0DhjduHEDp0+fRnp6OmxsbBAYGIj9+/cjIiKix7fkYBIs4vF4ADpGPCmKUow+maxhs1gsPH/+vK3dIhC6JIx2GBUUFOCvv/7CX3/9hbKyMowYMQI//PADIiMjFetjN2/exKFDh9rU2c4OkyLI8pFne653vsjIkSPx008/oU4ggFAoVGsbHByM8+fP44033iB93QmEl2AknrGxsQgICMD8+fMxfvx4mJubt7AJCQlBSEhIqzvYlWASSJOPPLURz/Lycty4cQONjY3Q19eHiYkJBg8erJOgGRsbIzY2Ftez7qMG1kCAaltra2vcvHkTdXV1ard1Egg9EUbieeLECfTp06etfenySCQStZH2yspKVFZWwtjYGHZ2doyumZGRgdraWkRHR8PZ2Rn6+vp4+PAhjh8/jtdeew0cDkdrPydMmIBjC4+h8loR6kYtBJfLVWnL4XBQWFjYrmlVBEJXgNGa58SJE1FUVNTi/ZSUFERHR7e6U10VsVisNtL+4nonk1Hj8+fPUVVVhVGjRiEwMBDW1tYwMzPDoEGDMHfuXCQnJ+u0jdLd3R1BQUEQCoVITExUaxsaGorz58+T9DQC4SUYiaeq7YYWFhYkoPD/kSfIMxVPJly7dg2DBw+Go6NjixGtg4MDZs6cifPnzzNuOPciLrP9EfBuOI4fP652OymHw4FIJCLZFATCS6idtm/fvh1AU5T2119/hampqeKYRCLBuXPnEBgY2LYedhFqa2thYGCgMccTAKO2zTweD66urujTp4/Ka7q6uuK1117DhQsXMGzYMI3V61/E3Nwc+vr6yMvLw927d9G/f3+VtjY2Nnjw4AGGDBnC+PoEQndHrXjK+xPRNI07d+40S+o2NDREREQEFixY0LYedhH4fD4sLS3VrnnKxdPV1VXj9Z48eYLRo0drbHPi7++P2tpa3Lp1C6GhoYyDSBRFwcq6KSf12LFjasWzX79+SEpKQkhISLsn9hMInRW14rl3714ATSPQuXPnNht5EprD5/NhY2Oj1qagoABAU2V+ddTX18PQ0JBxUCksLAyPHz9GaWkp43MAwNrKGmw2G6mpqSgrK1NZAFk+on3+/Dns7e2V2hAIPQ1Ga57Lli0jwqkBoVCodlRWXV2NqqoqGBkZaazSfu/ePQQEBK/ROToAACAASURBVDBOD6IoCq+//jpu376tVQBJX19fsd89Pj5erW1gYCAuXryosdwegdBTUDvy7Nu3L5KSkjQmSSclJbW6Y10NTUEb+ajT1dVV49RaIBDA3t5eYx+kFzEwMMDYsWORkpKCwYMHazzX37wp9Sx0kjcuXLiAkydPYvbs2Sr/ADg4OODevXsQCATkDymBAA3i+cUXX8DCwgLLly9vL3+6LDKZTO2ojOmUXSAQwNDQENbW1lr74OPjg2vXrqGyslLj+YNsggAAtDUNT09PPH36FJcvX1Y0jFOGhYUFHj58iIEDB2rtG4HQ3VArnvKGb6Txm2Y0JcjLg0WaxDMrKwt+fn467+iZNm0afv75ZwwbNgz6+voq7cSNEgCAAUsfkyZNwvfff49jx46pFc/Q0FAkJiYiKCgIhoaGOvlHIHQXVIrn9evXGV+EjESaxFNdqtCL03Z1CAQC2NnZaZV29CJcLhe9e/dGcXGx2nv9zmuqKP+W1zRER0cjLi4OWVlZePz4scrdZCwWCywWC6WlpRr/CBAI3R2V4jlnzhxGF6AoStFWoicjlUrVCp58h5azs7Pa6zQ2NsLCwuKVfHnttdewY8cO2NvbM0ot4nA4GD16NA4fPoxjx47hH//4h0rbsLAwnDlzBvPnz9dqTZZA6G6oFE95kzKCZgQCAfT09FSKp0wmQ3FxMYCmwIsqnj17BhsbG7V7zZnAZrMRGhqKnJwc+Pn5MTpnwoQJOHz4MJKTk7Fo0SKlxV8AwNTUFAKBAAKBoF06fxIInRUydGgFCgsLYWFhoXIkVlZWBolEAmtra7WFPB49egRvb28YGRm9sk9DhgzBs2fPUFdXx8je2dkZAwcOhFgs1pi21KdPH6SkpLyyjwRCV0ZtwGjdunVYs2YN9u3bp/Yiy5Yta1WnuhqaEuQLCwsBAI6OjmqvIxKJYGdn1yq1MymKwvjx45GYmIhBgwYxWkOdMWMGrl+/jiNHjmDatGkqRdzDwwPnzp2DSCQigSNCj0WtePL5fMhkMsU2TWWQIrlAXV2d2txH+XqnOvFsbGwEgFbNoXR1dYWBgQEqKytbiHt/y5Ztj4ODg+Hj44Ps7GwkJCSozLKgKArm5ubIycmBv79/q/lLIHQlGG3PlH8lKEcufKqQi6eTk5NKm8LCQjg6Omrcy64tU6ZMwe7duzF8+PBmqUvBVi3XQimKwt/+9jd88sknOHToECZMmKCyrUhISAiSk5Ph6+urc2YAgdCV0WrNMzMzE+fOncPZs2cVU1GC5t1FTKbt+fn5cHV1VZubqQtcLhd9+/ZVpErJqZMKUSdt2YYjPDwcrq6uKC0tRXJyssrrykVVUw94AqG7wkg8s7KyMHbsWMyYMQOffPIJPvroI8TGxuL//u//IBAI2trHTo9UKlW7fMFk5CmVSlVGuF+VkSNHgsfjQSwWK947mHcKB/NOtbBlsViYNWsWAOCPP/5QO6oeNGgQEhISyH53Qo+EkXiuXbsWffv2RWpqKlJTU3H9+nWcPHkSfD4fGzZsaGsfOz1SqVRlpL2xsVEhnurSlKRSaZvtGacoCiEhIXjy5Akj++joaPTq1Qv5+fm4cuWKSjszMzNUV1eTQsmEHgkj8czNzcWiRYtgaWmpeM/DwwNr1qxRO7XrKahrv1FRUQGxWAwLCwuVWy5pmgZN062+3vkigwcPRklJicaOmUDTlHzGjBkAgAMHDqgdWXp7e+PixYut5ieB0FVgJJ7+/v4t1syApvW0nl4ct76+Hmw2W6V4MlnvrKiogLW1dZum/VAUhREjRuD+/fuMptljxoyBhYUFsrOzcfv2bZV2Hh4eyMnJYZxPSiB0FxjtbR8xYgS2bdsGS0vLZh+8pKQk9O7du2097OQUFRXB3NxcpXgySVN6+vQpvL292zxnsm/fvrh8+TKjdWojIyNMmTIFv/zyC/744w8MGDBAqR1FUXBwcMDNmzcRFRXV2i4TCJ0Wrfa2v/HGGy3e6+l5nnw+H1ZWVip/D/KRp7pgkVAohKWlZbv8LidNmoRjx44hJCAALA33mzhxIv744w/cunUL2dnZKtsPBwUFISEhAaGhoW269EAgdCbI3vZXRCAQoFevXiqPMxl5ymSydtsnbmdnB2NjYzjTtrCytFJra2pqigkTJuDgwYM4cOAAPvnkE6V2FEXB1dUVaWlpakvaEQjdCbK3/RXR1M+cSZoSTdNq97y3NpMnT0banWuoEFZqtJ02bRr09fWRmpqqqEmqDH9/f9y5cwf19fWt6SqB0GlRu8PoRY4ePYrbt29DIpG0OLZ58+ZWdaoroU48aZrWGDCqrKyEmZlZqxQDYYqJiQlybctRwEvAYr+/q7W1trbGqFGjEB8fj//+979YtWqVUjuKouDh4YHLly9j1KhRbeE2gdCpYDTy/P777/HZZ5+hqqoKJ06cUOT1JSUlobq6uk0d7OyoSyKvrKxEQ0MDTE1NVU7LHz9+DC8vr3YVT6BJFOvq65slzqti5syZYLFYSExMxLNnz1Ta+fr64v79+2T0SegRMBLPY8eO4ZtvvsHWrVthamqKd999F5s3b8aSJUs0doJ8kdLSUixcuBBBQUGIjIxUW61Jk21iYiJiY2MRGBiIGTNm4NGjR82Oi8ViHD9+vM0jwBKJRGOwSN16p1AohLW1dbsH3iiKahqB5uZqtHV0dMSwYcMgk8nw559/qr2mt7c3yf0l9AgYiWdFRYUiJcnCwgLl5eUAmmpGnjt3jvHN1q9fD0dHR6SmpuLHH3/E9u3bVeYQqrMtKSnB6tWrsX79ely9ehWxsbFYuXKlIo2qsrISwcHBWLNmDWPfdEVdBfnOFix6GQtzc/D5fDQ0NGi0/dvf/gYAOHXqFKqqqlTaeXl54dGjR2T0Sej2MBJPCwsLhRD0799fIZjl5eUai2LIqa2tRUpKCj744AOYmpoiMDAQ48aNQ0JCgta2iYmJCAsLw7Bhw8DlcrFgwQKUl5cjOzsbAGBpaYnMzEzs2bOHkW+6QtO02sZvTNKUAHRoek9UVBQyMzM1Js57enoiLCwMIpEIR44cUWlHURT69u2Ls2fPtrarBEKngpF4RkZGIjMzE0BT9PXo0aOYOnUqVq5ciWnTpjG6UV5eHoyNjZv153F3d1c6bdRky+PxmvUCYrPZcHFxAY/HY+RLa1FTUwNDQ0ONI09V4ikWi6Gnp9chBYXH+4zEeJ+RCAwMRH19PaP96bNnzwbQFDxUN/p0c3PD06dPSdEYQreGUbT9iy++UHwfGhqK/fv3Iz09He7u7hg5ciSjGzU0NLQIihgZGSmd3mmyFQqFLYr7qrqWJrZt24bt27drfR7QtC5ramqqcuSpqSAIn8+Hg4NDh2xxDXUKUnw/bdo0HDx4EFFRUWqbuvn5+WHQoEFIT0/Hf//7X7zzzjtK7SiKQlBQEE6ePImZM2f2+I0UhO6J1nme1dXVcHV1xcKFCxkLJ9Akbi9P8RsaGpTmN2qy5XA4LVKmVF1LE8uXL0d2dnazV1JSEqNzS0pKYG5urlJw5E3fVK15FhcXw8XFpUOKCRfVlKCopgRAU+TdzMwMz58/13jevHnzADQFEeVr38pwcnICn89Xa0MgdGUYiWdNTQ02bNiAkJAQDB48GEOGDEF0dDQOHjzI+EZubm6oqalp9mHi8Xhwc3PT2tbd3R05OTmKYzKZDPn5+XB3d2fsT2tQW1urcr2ytrYWtbW1MDIyalaN6kUkEkmb1fDURNyN/Yi7sV/x89SpU3H79m2NSf/e3t6IiIiAWCzG/v371doOHToUf/75p8ZrEghdEUbiuWrVKty6dQtffvkl4uPjcfToUbzzzjvYunUr4ymvqakphgwZgq1bt6K2thb37t1DfHw8Ro8eDQCYO3euIh1Jk210dDSuX7+Oixcvoq6uDnFxcbC2tkbfvi378rQlMplMZaBFPuq0t7dXOW2VyWQqy9S1NxwOB15eXow6BMybNw8URSlquqrC3NwcFEU1+0NHIHQXGIlnWloaNm7ciNjYWHh5ecHX1xevv/46NmzYoLGz5ots3LgRBQUFCA8Px5IlS7B06VKEhoYCAAoKClBZWcnI1snJCZs3b8bnn3+OsLAwJCUlYcuWLQqRqqqqQlhYGN59912UlZUhLCwMq1evZuwnU9SNqDRN2eW0d3K8OkaPHo3s7GyNifPu7u4YNWoUZDIZfvrpJ7W2Q4cORXx8vNKdaQRCV4ZRwMjBwUFpLqCnpyejHSovXueXX35ReuzlxGp1tkBTvckxY8YoPWZhYaG242drwUQ8VQWL5PmhnakeKpvNxqBBg/D48WONXTHnzZuHCxcuICUlBXfv3kW/fv2U2unp6cHNzQ3nz5/Ha6+91hZuEwgdAqOR55IlS5R20Lx//z4CAgJa3amugjrx1BRpLyws7LBIuzrCwsLw7NkzjZkLNjY2eP311wEA//nPf9RuU/Xz88Pdu3dRU1PTqr4SCB2JypGnr69vs7U6mqZbrCnSNK2yNW1PQF3vIk3T9sLCQgwYMEBtalBbMtVP+aidoihMnDgRJ0+eRGRkpFr/pk+fjvj4eDx69AhJSUmIjY1Vec2wsDAcPnwYc+fO7bBnJhBaE5XK99tvv7WnH10SiUSiMs1I07RdIpG0WcM3JgTZqw6uubi4wMrKCsXFxWp3R3E4HCxYsABff/01du3ahfDwcHC5XKW2tra2uHfvHnJycuDt7f3K/hMIHY1K8Rw0aFB7+tHlkEgkoGla6ShKKpWitLQUFEXBzs5O6fkymUyl0LQHuZVNPancLV2UHp82bRq2bdsGGxsbtTugRo4ciRMnTuDBgwfYu3cvFi9erNI2MjISf/31F5YtW0YqzhO6PIznT1lZWVi1ahWmTp2KqVOnYtWqVcjKympL3zo1FRUV4HK5SkeepaWlaGxshI2Njco1TZqmOzTSvuf2Iey5fUjlcTabjZiYGI0N41gsFlasWAGKonDkyBG1VZr09fURHBxMcj8J3QJG4pmQkIDXX38dLBYLU6ZMwZQpU8BisTBr1iycOnWqrX3slJSVlancmqmpmpJMJgOLxeqQPe3aEBgYCKFQqHHfe58+fTB+/Hg0NjZi69ataoNHzs7OEIlEuHv3bmu7SyC0K4yiPVu2bMHatWsVhSHkBAUFYcuWLRg7dmybONeZke9rV5YAr2m9s7i4GHZ2dp1ePIGmlh2HDh3CsGHD1AZ63nrrLVy6dAkZGRk4efIkJkyYoNI2PDwc8fHx8PDwULn7ikDo7DAaeRYWFipdAx04cKBilNXTaGhogL6+vtJjmsSzoKAArq6uHbKnXVtsbW1ha2uLkpIStXZmZmZYsWIFACAuLk5txXkWi4Vhw4Zh7969JHme0GVhJJ5eXl5Ki2WcPXsWnp6ere5UV4BJjqeqabtYLO6wAsi6MHnyZGRkZGgUumHDhiEyMhJCoRDffvut2rVSc3NzuLq64vjx4xpriRIInRFG0/bVq1djyZIlSEtLg5+fH4CmBPmMjAzs3LmzTR3srNA0rXLPuqaRZ0dH2gHgb0GTGNvq6+sjOjoad+/eRUhIiNoScytWrMDdu3dx69YtnDp1CuPGjVNp6+3tjatXryItLQ3h4eFa+U8gdDSMRp7h4eE4c+YM+vfvj4KCAhQUFCA4OBinT5/usf/ppVKp0hETTdMaxbOxsbHDU3V8bLzgY+PF2D4oKAgNDQ0adwlZWlpi+fLlAJp2HqmbvsuT569evaq2rTGB0BlhJJ5r166FRCLB+++/j+3bt2P79u344IMPVIpDT0BV+5GqqirU19fDxMRE6dRcHmnv6G2Z2eVPkF3+hLE9RVF4/fXXcfXqVY1pRsOHD0dERATq6+vx/fffa0x1iomJwf79+9VWpycQOhuMxDMpKYmsS72Eqt1FL27LVBWJt7e37/BI+4GMYziQcUyrc8zMzODv74/Hjx+rtaMoCitWrICpqSlu3LihtE/VixgYGGD48OHYvXs3o2Z0BEJngJF4DhkyBDdu3GhrX7oUqva1a2o3zOfz4ezsrDJS39mJiYlBfn6+RpGzsrLCsmXLAAA//vijxhbH5ubm6N+/P3799VfGTQUJhI6EUcCod+/e+OqrrwBAqWBMnjy5db3q5MiFQ9nIk0mkvSP3tL8qLBYL06dPx7FjxzTmfkZHR+PGjRs4d+4cPvvsM+zYsUNtqxRHR0cIBAIcOnRIsSmDQOisMBLPv/76C1wuV2nVeIqiepx4lpWVwczMTKfdRVKptNNUj9cVJycn2NnZoaCgQGkbFTny6fujR4+Ql5eHLVu2YM2aNWqj9d7e3rhz5w5OnjyJ8ePHk+ZxhE4LI/F8uVBxT+fZs2cqt2ZqmrbTNN3haUqtwaRJk7Bt2zb06tVL7WiSw+Hg448/xtKlS5GYmIh+/fqpLGItp1+/frhx4wbOnj2L1157jQgooVOicV6Ul5eHPXv24NChQ2rTTnoS8pGnuq2ZysRT3qe9oyPtADAveAbmBc/Q+Xw9PT3MmjULV65c0Rh9d3Nzw8qVKwE0tXp++vSpWnuKohAaGori4mJcvHhRZx8JhLZErXjevHkTU6ZMwZYtW7B582aMGTOGBI6gOlgkEAhQXV0NQ0NDWFtbtzguDxZ1dKQdaCpFp6ocHVMcHBzg4+OD7OxsjbaxsbEYM2YMxGIxPv30U9TV1am1pygKAwcORE5ODi5fvvxKfhIIbYFa8dy6dStiYmJw48YNXL9+HRMnTsQ///nP9vKt06KqatCLrTdUjUrt7Ow6RaQ9o+QBMkoevPJ1Ro4ciZKSElRXV2u0XbZsGTw8PFBYWIgvv/xS44iVxWJhyJAhyMrKQlpa2iv7SiC0JmrFMyMjA3PnzgWbzQabzcYHH3wAHo+HsrKy9vKvU6Kq5bBcPFVVX5dKpSqn++3NkawEHMlSn3/JBIqiMGfOHFy9elVjipGhoSE2bNgAU1NTpKWlYdeuXRqvz2KxEBERgTt37pARKKFToVY8hUIhrKysFD+bmZlBX19f45Sru6NKJJjU8ezKaUqqMDMzQ1hYGO7du6dxM4WzszM2bNgANpuNP//8EydOnNB4fRaLhaioKDx9+pQUEiF0GjRG23/99ddmH/jGxkbs27cPFhYWivfkydA9BYlEolOaEk3THb6nva0YPHgwsrOzUVJSonHbbv/+/fHee+/h22+/xbZt22BhYYHIyEi157BYLISGhuLRo0fYs2cP3nzzzS5R0o/QfVErngMHDmzRamPAgAHNAgSdYQra3sh7rr+MujQlgUAALpfbKYJFbQFFUZg9eza2bt0KMzMzjelYY8aMQXl5OX799Vd88cUX+OKLLxAcHKz2HBaLBR8fHxQVFWH79u2YP39+txzJE7oGasVTWa/2no5IJAJN00rFs6Cgqamas7Nzi2M5OTlwd3fvtuIJNJWue/PNN/H7778jJiZGY1vqN954A9XV1fjrr7/w8ccf45tvvoGvr6/acyiKgpOTE7hcLnbt2oXp06fD1dW1NR+DQGAE2f+mJap6F9XW1qKqqgpGRkawtbVtcV5NTQ169erVabYcvh36d7wd+vdWv66trS0iIyNx48YNjWuTFEXh3XffRXR0NIRCIdasWYOHDx8yuo+FhQVGjBiBEydOID09vTVcJxC0onN8krsQJSUlSrdmykedTk5OKtsRd6YppqOZPRzN7Nvk2iEhIbC0tGQkhCwWC6tXr0ZkZCTq6uqwevVqxl1ZjYyMEBUVhQcPHpCOnIR2h4inlpSXlytNN+Lz+QAAFxflieeNjY2dak/7jcIM3CjMaLPrT5w4Ec+fP2fU40pPTw8fffQRhg8fjvr6eqxZswZ37txhdB89PT0MHDhQUXuhsrLyVV0nEBhBxFNLVEXa5ZXQlYmnWCyGgYFBp1rvjM9ORHx2Yptdn8Vi4c0338Tdu3cZJdDr6elh3bp1iImJgVAoxLp163D+/HlG96IoCu7u7ggPD8eePXvILjhCu0DEU0tU7S5SN/LMzc2Fs7MzjIyM2tS3zoaBgQEWLVqES5cuob6+XqM9m83GmjVrMGXKFEgkEmzatAmHDx9mfD8TExOMHDkSWVlZ2L17d4/PRya0LUQ8tUTVupq6keezZ8/g7OzcI/MSTU1NsWDBAiQmJjJqM8xisfDuu+9i0aJFAICdO3fiP//5j8o/Wi+jp6eH0NBQ+Pn54eeff0ZycjJJqie0CUQ8tUTZ1kyZTKZY21OWptTZgkXtjbW1NWbNmoXTp08zCurI+yWtWbMGbDYbhw4dwmeffQaBQMDofhRFwcLCAsOGDUN1dTW2b99OGswRWh0inloiFotbjCBLSkoglUpha2urtLalTCbrVMGijsDNzQ1jx47FuXPnGI8iY2NjsWnTJhgbGyMlJQVLlizR2D/pRfT19eHr64vw8HCcPXsW+/btg1Ao1PURCIRmEPHUApqmlTZ+y8vLAwClydoikQiGhoZqCwZ3BMsGz8OywfPa9Z4BAQGIiorChQsXGAtoaGgodu7cCS8vLxQXF2PFihU4ceKEVlNxY2NjDBkyBG5ubvjpp59w5MgRRksIBII6iHhqQWVlJTgcTotou7y4r6enZ4tznjx5And3904XLLIxtoKNsZVmw1ZmwIAB6N+/P1JTUxkLoJOTE7Zt24bx48dDIpFgy5Yt+OKLLxgFoeSwWCzY2Nhg+PDhsLKywq5du3D8+HHSbI6gM0Q8taCoqAhmZmYtRp5y8fTw8GhxTkVFBWxtbTtFDc8XuZJ/A1fy2z+lh6IohIeHo0+fPjh//jzjEaiBgQHee+89fPjhhzAyMsL58+exePFi3L59W6v76+npwc7ODkOHDoWJiQni4uJw8uRJkmBP0BoinlpQVFQECwsLrUae8hqenY2zOZdwNudSh9yboigMHz4cQ4YMwalTp7SaQkdHR+PHH3+Ep6cnioqKsGrVKnzzzTeoqanRygcDAwM4ODggPDwcHA4H//73vxEfH0+m8wTGEPHUAnkPohdpaGhAYWEh2Gy20jXPxsbGTimeHQ1FURgwYACmTZuG+Ph4xpF0oGlteceOHXjrrbegr6+PM2fOYP78+Th//rzWaUmGhoZwcHBAREQEOBwO4uLisH//fkWFLAJBFUQ8tUDZ+lhubi5omoaLi0uLxm51dXXgcDjdtoZna+Dp6YlFixbh/PnzKCkpYXyevr4+Zs+ejbi4OAQFBaGqqgqbNm3CRx99pAjgaYOBgQEcHR0RFRUFLy8vJCcn45dfftFqaYHQsyDiqQVSqbTFyEbdlP3x48fw9PTsdMGizoatrS2WLl2Khw8fMmom9yIuLi7417/+hffffx9cLhfp6elYtGgRvv32W5SXl2vti56eHszMzDBgwAAEBgZCIpHgp59+wsGDB1FaWqr19QjdFyKeWiCRSFoUBFEnnlVVVejVq1eP3FmkLVwuFwsWLIBIJEJKSopWoz0Wi4Vx48Zhz549mDhxIgAgISEBc+fOxe7du1FVVaW1PywWC1wuF87Ozhg4cCBcXV2RkJCA3377DSkpKWTXEoGIpzYoqyCvTjxlMhnMzc3bxTdt+WDo2/hg6Nsd7UYzDAwMMGPGDISEhCA+Ph7Pnz/X6nxLS0usWLECu3fvRkREBEQiEQ4cOIDZs2dj27ZtWi0LyKEoCkZGRjA3N0dISAj69u0LgUCAf//73zhw4ADS09PJtL6HorGHEaGJ6upq6OnpNQsYNTY2IicnB0BL8ZTJZKAoqtNuyzQz7Jw7nlgsFgYMGAAvLy/8/vvvMDc3R3BwsFbtXlxcXPDJJ58gKysL+/fvx9WrV3Hs2DGcOHECI0aMwMyZM+Hl5aW1b2w2GyYmJjAxMYGDgwPEYjHKysoQFxcHU1NT9OrVC5GRkWSZpodA0WT+0QI+n4+YmBgkJSUp9qrfv38ffD4fLi4uilQlHo+HRYsWoVevXti/f3+zazx8+BAcDgcxMTGdLscTAC7wmvqgD/cY0sGeqEYikSA1NRXXr1/H4MGDYWNjo9N1eDwe/vvf/yI5OVkxSvT398eECRMQFRXVItCni58ikQgCgQCPHz8Gh8OBmZkZIiIiYG1t/UrXJnReyMiTIYWFhS1yPOUVz/38/FrYP3v2DEOHDu2Uwgl0DfHU19fH8OHDERQUhKNHj+LOnTs6iZ2HhwfWrl2LefPm4c8//8TZs2eRmZmJzMxM/Pjjj4iNjcVrr72m02hU7qe+vj5MTExga2sLoVAIiUSC48ePo7GxEVwuFxYWFggLC4OlpaVO9yB0Poh4MkRZQRB14imVSskHpZWwsrLC3LlzkZubi/j4eJiYmCAsLExjg7mXsbe3x7Jly7BgwQIkJyfjxIkTyMnJweHDh3H48GF4enpi+PDhGD58uMr20ZqQT+0BICwsDGKxGBKJBA0NDTh69ChomgaXywWXy8WAAQPg5OSk030IHQ8RT4YoS1PKzMwE0FI8BQKBIshAaB309PTQu3dvLFmyBNnZ2UhOToaenh7Cw8O1XmPkcDgYN24cxo4di+zsbJw9exbnz5/H06dP8fTpU/z888/w9vZGeHg4wsLC0Lt3b51abFMUBUNDQxgaGsLExATW1tYQi8UQi8WQSqW4cOEC6urqwOVyweFwYGpqigEDBpCpfheBiCdDXt629+zZM/D5fBgbG6N3797Njt2/fx9+fn4ae5cTtMfQ0BBBQUHw8fFBXl4ezp07B6lUioEDB2o90qcoCr6+vvD19cXixYtx48YNXLhwAWlpaXj06BEePXqEPXv2wNraGoMGDcLgwYMRGBio846xF8UUaOoAKpVKIRaL0djYCLFYjFOnTkEsFoPD4cDIyAjGxsbw9vaGu7t7p+m8SmiCiCcDaJqGSCRqNk2U98kJDg5uMX0UCASwt7cn/9nbEENDQ4WolJaW4syZM6ipqYGPj4/SAi2aMDAwQHh4OMLDwyESiXDjxg1cvXoV6enpqKioQEJCAhIS7dbFfQAAFy1JREFUEgA0bQ/19/dXvJydnXUemcrXS+VYWlpCIpFAIpGApmlIpVLcvHkTycnJCuE1MDAAh8OBj48PXF1dyf+zDoKIJwPKysrA4XCUimdoaGgz2+rqanC53E4/9VoXtayjXWgVDAwM4OzsjHnz5qGyshKpqak4e/YsWCwW/Pz8dFq7NDQ0xNChQzF06FDQNI2cnBxcu3YNN27cQHZ2NvLz85Gfn68QUzMzM/j7+8PPzw8BAQHw9vbWudkfi8VqNjoF/jdClb+AppnQ9evXkZiYCH19fcU5enp6MDU1hbu7OxwdHV85k4CgGpKqpISXU5UuX74MmqZhY2MDiqIgkUgwffp01NXVYe/evXBwcFCce/nyZfTv3x+hoaFkRNAB0DSN+vp6PH/+HGlpaSguLgaLxYKXl5dOI9KXEYvFyMnJQVZWFu7fv4/MzMwW7Y5ZLBacnJzg4eEBd3d3eHh4wMPDAw4ODq2620w+MpXJZM3W5EUiEcrLy1FTUwMWi6UY3crzlM3NzeHm5gYHBwetg26E/0F+cwyorKyEg4ODYmp27do11NXVKT4QchobGyEUCuHo6NjphfPM44sAgFF9hnWwJ60LRVGKaLazszPq6upQVVWFW7duITExERRFwcDAAH5+fjrNDuTn+vn5Yfr06aBpGiUlJQoxzcrKAo/HQ0FBAQoKCnDp0v/K/hkaGsLNzQ0uLi5wdHSEg4OD4mVlZaX1/xll036gqemelZWVQljlfbfk/3/r6uqQnp7eTFz19PSgr68PFosFFosFCwsL2Nvbw9bWFubm5jotS3R3iHgyQCQSNYu0nz17FgDw2muvNbO7c+cO/P39dU7mbk/SCm4C6H7i+SIURSl2BDk5OUEoFKKurg6lpaW4ffu2opCynp4ebG1t0bt3b62n2xRFKQQwJiYGQNPoND8/HzweDzweD7m5ueDxeCgrK1MEol7GwMAA9vb2cHBwgJ2dHaytrWFjYwNra2vF9yYmJoxFjMViqZyym5iYwMbGBo2NjQpxlW8eoCgKNE2jrq4O9+/fh0AgQENDg2LU+rLIstlssNlsmJqawtLSEhYWFjAzM4OJiUmnH0C8Ku0qnqWlpfjwww+Rnp4Oc3NzvPPOO3jjjTd0sk1MTMRXX32FkpIS+Pr6YtOmTfD29gbQlFa0adMmRa+bcePG4Z///KdOCevy3SPy6Rafz8fVq1fBZrMRHR3dzK6srEyn1BlC20NRFIyNjWFsbAxbW1v07dsXDQ0NqKurQ21tLXg8HtLT0yEWiwE0iY+xsTFcXV21Lu5iYGCA3r17t8jCEAgE4PF4KCwsRFFREYqLi1FcXIySkhJUVVUp1lJVYWhoCGtra1hZWcHc3BxmZmbNvsq/l//M5XJVCphc9FR9JkxNTWFnZ4fGxsZmr5e7x8rFVigUIi8vDw8fPkRDQ4PiM8Nms5uJrFx05S/5TMDCwgLm5uYwNTVVzBw6+5JCu3q3fv16ODo6IjU1Fbm5uVi0aBH8/f0RHByslW1JSQlWr16N77//HqGhofj999+xcuVKnDp1ChRFYc+ePcjIyMDJkycBACtWrMDPP/+Md955R2ufr1+/DhcXFxgaGoKmaezcuRONjY0YM2ZMs2nfhQsXMGjQIJ2Tqwnti1wc5WLq6emJqKgoiEQixVbLsrIy8Hg8ZGdnNyv+If/gW1pawtHREZaWloxGWSYmJggMDERgYGCLY/X19SguLkZRURHKy8tRXl6OiooKVFRUKL6vr69HUVGRos01k2eUC5GxsXGL71/8amRk1OxlaGio+CpPmzIwMFD5R0S+t18usjRNt/j+ZdGVI5VKUV1djaKiIohEIkgkEojFYtA03UJsXxRjiqJAUZRChFksluIcPT09cDgcRT1dIyMjxc/yZ3zVpYh2E8/a2lqkpKTgypUrMDU1RWBgIMaNG4eEhIQW4qnJNjExEWFhYRg2rGnKuWDBAuzatQvZ2dnw9fXFqVOnsGTJEtjZ2QEA5s+fj507d+oknkVFRXBxcQFN0/jpp59w7do1mJiYYO7cuQCaFu1TU1Ph6uoKHx8fnaOshI7nxW2W1tbWcHNzQ0hIiCJ1SCwWQyQSQSgUQigUorS0FLm5ucjIyFBaou7F0ZW+vj5MTU0VywhcLheGhoaKD7CxsTG8vLzUbhGtr69XCGpNTY3iVV1djerqasX38q/19fWora1FbW1tq/2ODAwMYGhoqPhdvfwyMDDQeEw+/ZePRlWNTuUC+bKAvmj7oj1FUWCz2YpAmkAggEwmg1gsVixPvPhis9ngcDiIjY3VadDTbuKZl5cHY2NjWFhYKN5zd3fH5cuXtbbl8XiKgh1A0xTExcUFPB4Pvr6+LY67u7sjNzeXsa/yZmCzZ89WpIaUl5ejrq4ObDYbb7/9Nurr63HhwgWUlZWhd+/ecHNzg1AoBJ/PZ3yfjkRQ0fSB6ir+dibkYmBmZgYvL69mEW/5V6lUqhBdmUwGoVAIgUCA4uJihfjKlwiYIh9pvYh8mu7i4gIAilGXXHhe9EMkEin+ALz8ku96ku+Aevklt5FfqzshX5rRdpmg3cSzoaGhxVqgkZGR0vaxmmyFQmGLoMyLx18+n8PhQCgUNos4ytm2bRu2b9+u1OeX6z/KdwzFxcWpfM6uxgns12xE6PEoi+p3FxYvXtysghpT2k08jYyMWvQAamhoAIfD0dqWw+G0+Ov34vGXzxcKhSrXOJYvX47ly5e3uFa/fv1w9uzZHlcFXp7f2hMhz95zn93e3l7r89pNPN3c3FBTU4Py8nLFqPH/tXfvUT3ffwDHn/2asoakI7bOSJpvRRe3bhPVOigppw2JXFJzaXMkdtTIbVhSh2JroTG3c2iFxYjQISoUhXVDaLO0hnTZ16rv7w/H9/gqU1/lu9X78Vd9vu/P5/N696lX78/t9b516xa9e/dudlsDAwOF5+fq6uq4c+cOBgYGwNMSZIWFhXzwwQf/uJ+XeTZqbc46bUlz/wO3JaLv7ZMyd/bf2INYnTt3xtbWlqioKB4/fkxubi5JSUmMHj0agGnTprFr164mtXVycuLChQukpqZSVVVFbGwsurq6mJiYADBq1Ci2bt3K77//TmlpKdu2bZOvKwiC0BLe6KNKX331FSEhIdjZ2aGtrU1AQID83fC7d+8qvOb2T2319fVZu3Ytq1atkj/nuXHjRvlpua+vL7/++itjxoxBJpPh4uKCn5/fm+yqIAhtnHi3/SUkEkmzp8FtC9prv0H0XfS9edSXL1++vOXDaRusra1VHYJKtNd+g+h7e6VM38XIUxAEQQlt+819QRCEViKSpyAIghJE8nzB/fv38fPzw9zcHHt7e/njU21ddHQ0xsbG8lqVpqameHp6qjqsVlNZWcmqVatYuHChwvJLly4xduxYBgwYgJubG5mZmSqKsPW8rO8+Pj6YmJgo/A58/fXXKoqyZcXExODi4oKlpSWOjo5s27ZN4XOljrtMUODv7y9bunSprKKiQpaTkyOztraWZWVlqTqsVhcVFSULCgpSdRhvxLFjx2QmJiYyiUSi0Ofq6mqZjY2NLD4+XlZVVSVLTEyUWVlZySorK1UYbct6Wd9lMplsypQpsn379qkostYVExMju3Llikwqlcry8vJkdnZ2stOnT8tkMuWPuxh5PudZNacFCxY0qOYktB0jR47k+vXrBAQEKCxPT0+ne/fufPzxx2hpaTFu3Dj09PQ4f/68iiJteS/re1s3a9YszM3N0dDQQCKRMGTIEPLy8gDlj7tIns95WTWn5lRk+i87cuQIAwYMYPjw4Xz11VfNrvrzX/diNS5oX8cfYPny5VhYWDB69Gj27t2r6nBaRV1dHVevXpWX/1P2uP+7SzW/Yc2p/NTWeHt7M3XqVDp16kRRURHBwcFEREQQHBys6tDemJqamgaFatrL8QcICwtDR0cHdXV1MjIyWLBgAXp6evLpRdqK9evX061bNxwcHADlj7sYeT6nOZWf2hpdXV20tbVRV1dHIpHg5+enUHylPXhVta627r333uPtt99GQ0MDe3t7XF1d29zvQGxsLCdOnGDz5s3yYiDKHneRPJ/zfDWnZ5pbkamtqK6u/k9MZNeSDAwMKCoqUlh269YtebWu9qaxurn/ZZGRkRw8eJBdu3ahp6cnX67scRfJ8zmvqubUlq1bt46cnBykUil5eXls2bIFd3d3VYf1RllbW1NeXk58fDzV1dUkJiZy//59bG1tVR1aq7t//z6RkZHcunWLJ0+ecPLkSU6cOIGrq6uqQ2sRQUFBZGdns2fPHvn0PM8oe9zF65kvuHfvHiEhIVy8eBFtbW38/f3l8xW1ZWvWrOH48eOUl5fTs2dPJk2axPTp09vkfN0nT54kODiYv/76i/r6erS0tAgMDMTLy4vMzExWrlxJcXExvXr1IjQ0FBsbG1WH3GJe1ndXV1cWLVpEbm4u1dXVGBoaEhgYiL29vapDbhESiaTRwubXr18HUOq4i+QpCIKgBHHaLgiCoASRPAVBEJQgkqcgCIISRPIUBEFQgkiegiAIShDJUxAEQQkieQoq5ePjQ3R0tKrDaHElJSVIJBJKSkr+sV1tbS2enp4Namsqs63XJZPJmD17NnPmzGnV/bQVInm2U5WVlaxfvx5nZ2fMzMwYMWIEM2bMICEh4bW2m5CQgEQiaaEonxZpdnJyalZ7Hx+fFtt/a9u+fTuVlZWsXr1a1aGgpqZGeHg4ly9f5siRI6oO519PVFVqp+bMmcMff/xBSEgIffr0oby8nDNnzrBv377XqiA/atQorKysWizOadOmMX78+Bbb3r+JVColLi6O0NBQNDU1VR0O8PQV5c8++4yYmJg282pmaxEjz3aovLyczMxMgoKCcHJyok+fPgwZMoTAwEDi4uLk7WpqaggNDcXa2hpra2sWLFjAgwcP5J9LJBK2b9/OzJkzGThwIBs3buTYsWNMnTpV3mb37t189NFHmJmZYWNjQ2BgIH/++WeTY92xYweLFi2Sf+/j48OSJUv44osvsLW1xdHRkeTkZAAyMjLYtGkTmZmZSCQSJBIJGRkZAKSlpTFu3DjMzc1xdXVVKHAdHR2Nl5cX33zzDc7OzgwYMIAdO3Y0GPE+fvwYMzMzLl++DICvry+2traYmZnh7OzMjh07mtwvgJSUFNTU1Bg5cqTC8traWiIjIxk6dChWVla8ODt4VlYWbm5uDBw4kEGDBjF58mSuXr0KQEFBARKJhNu3byus4+npyZYtWwA4f/4848ePl09JERwcTE1Njbyth4cHxcXF5OTkNKs/7Y1Inu1Qly5d6NixI+fOnaOurk7hMy0tLfnXy5Yto6SkhNjYWLZu3UpZWRkrVqxQaL99+3YmTJhAYmJioyPEzp07ExQURGJiIps2baKwsJCwsLDXiv/w4cOYmpqyc+dORo8eTUhICNXV1VhaWjJ16lQsLCxISUkhJSUFS0tLbty4weeff463tzcJCQnMmDGDRYsWUVBQIN9mdnY29+7dY9OmTcTHx+Pm5kZpaSlZWVnyNsnJybz77rtYWloCYGZmxoYNGzh48CCzZs0iLCyM9PT0JvcjPT0dZ2dn/vc/xT/DDRs2cPjwYSIiIti9e3eDd6zV1dWZNGkSe/bsYc+ePejq6jJv3jzq6+vp168f/fv358CBA/L2hYWF5OXl4eHhwcOHD5k7dy6urq78+OOPLF26lIcPH1JVVSVv36lTJ+zs7JrVl/ZInLa3Qx06dCAkJISVK1dy5MgRzM3NMTU1ZdiwYQwZMgR4epPi6NGjpKWl0blzZwAWLlyIt7c3dXV18iILa9eu/cfqMy9WZpo4cSI7d+58rfhnzJjB9OnTAfD39ycuLo7i4mJMTU3p0qULmpqaCpXBt27dire3NxMmTADAyMiIY8eOkZycTL9+/QAYPHgwq1atUtjPhx9+SFJSEoMGDQKeJu3n+xMYGCj/2tDQkL1795Kdnd3kQiJFRUW4ubkpLJNKpfzwww9ER0czfPhw4Gm9yfDwcHkbCwsLLCws5N8HBATg7u5OWVkZPXr0wNPTk7i4OObNm4eamhoJCQnY29ujp6dHTk4O1dXVuLi40LNnT/r27dvoNWUjI6MGZdoERSJ5tlMTJ07EwcGBU6dOcf36dVJTU/n2229xd3cnPDycwsJCpFIpdnZ28nVkMhm1tbWUlZXRs2dPgAajpheVlJSwc+dOMjMzKSsr49GjR3Tv3v21Yn++0pOOjg7w9JT6ZfLz88nPz1c4ra6trVUoTdZYxR13d3dWr15NSEgIjx49IiMjQ2HknZKSwuHDh8nNzaWiooKKigocHR2b3I+HDx82KI92584dpFKpPGE3pqqqil27dnH27Flu374t7/uzgr5ubm6EhYWRkZHB0KFD+emnnwgNDQXA2NgYY2NjPDw8sLe3Z+DAgbi4uNCtWzeFffTo0YP8/Pwm96U9EsmzHevRowdeXl7y7xMTE1m8eDGzZ8+mrq4OLS0t4uPjG6ynq6vbpO3X1NTg5eVF7969+fTTT9HX1yctLY39+/e3WB+aUjKvvr4eX19fxo0bp7D82Yj6ZZydnQkNDeXcuXPcvXsXc3Nz3n//fQCOHz9OYGAgU6ZMYdmyZejo6LBs2bJmxS6TyRokbalUCiCvct6YhQsXkp+fz6xZszAyMqK6uho/Pz/55127dsXR0ZHExESkUim1tbXypK6hocH+/fs5ceIE6enpfP/992zevJlDhw4pFD6WyWRtshxhSxLJsx169OgR9fX18lHbM0ZGRsDTEcyzP8onT55gYmKi1H4KCwspKyvj6NGjdOrUCaDVTwU1NTUbTKnQt29fbty4IZ/wq6k6duzIyJEjSUpK4u7du3h4eMg/S0tLw8HBgcWLF8uXNXe6Dh0dHUpLSxWW6evro6amRlFREWZmZo2ud+7cOSIjI+VzCzX2/KenpyeBgYE8ePAAd3d3OnToADwdcWtoaODq6oqrqytSqRQbGxsuXryoUPS7tLRUYSJEoSFxw6gd+u2333B1dSU6OpqsrCxu3rzJqVOnWLJkCcOGDcPY2BgDAwNGjRrF/PnzOXnyJDdv3uTEiRPMnTu3yfvR19enQ4cOxMfHc+PGDZKSkoiNjW3Fnj299njt2jVSU1PJzs6mrKwMf39/Tp8+zbp16/jll1+4du2a/MbQq3h4eJCcnMy1a9dwcXGRLzcwMCA3N5eLFy9y9epVoqKiuHLlSrNibey6oo6ODiNGjCAsLIzc3FzS0tIaTMJnYGDA0aNHKSoq4uzZsw3uxgPY29vzzjvvkJqaqvDo2aVLl/Dz8+Ps2bMUFxfz888/8/fff8uv/T5TUFAg/2cqNE4kz3aoV69e+Pr6cubMGQICAhg7diwrVqzAwcGBqKgoebs1a9ZgZ2fH0qVL8fDwYP369RgbGzd5P7q6uqxevZq4uDg++eQTDhw40KxrgspwcnJizJgxzJ8/n9mzZ1NaWoqxsTExMTFcuHCBCRMmMHPmTK5cuYKpqekrt2dtbU2XLl0YPnw42tra8uXe3t5YWVnh7+/PnDlzkEqlzX45wMbGhuPHj/NiPfLly5ejpqbG5MmTCQ8Pb/AzW7NmDQUFBXh6ehIRESG/sfQ8dXV1nJ2d6d+/v8IxMzQ0pGvXrnz55Zd4eHgQFxdHZGQkhoaG8jYVFRWcP3++TVXQbw2ikrwgqIhUKsXR0ZG1a9cyYsSIFt/+2LFj8fLyYvLkyc1ab9u2bRw8eJBDhw61eExtiRh5CoKKaGpqMnPmTCIiIho8b/u6Ll26RHFxcYNHoV7lwYMHfPfdd+L99iYQyVMQVGjatGm89dZbLf5u+759+3B2dla41PAqdXV1zJ8/n8GDBytc3xUaJ07bBUEQlCBGnoIgCEoQyVMQBEEJInkKgiAoQSRPQRAEJYjkKQiCoASRPAVBEJTwf4N+GyjPeHEjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 324x252 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(4.5,3.5); \n",
    "\n",
    "clrs = ['k', 'C7', 'C1']\n",
    "\n",
    "disease = 'COVID-19'\n",
    "Ymedian = [np.median(yy[j]) for j in range(len(x)-1)]\n",
    "Ylower = [np.percentile(yy[j],[2.5])[0] for j in range(len(x))]\n",
    "Yupper = [np.percentile(yy[j],[97.5])[0] for j in range(len(x))]\n",
    "ax.plot(x[1:], Ymedian, lw=2, color=clrs[0], label=disease)\n",
    "ax.plot(x, Yupper, lw=.3, color=clrs[0])\n",
    "ax.plot(x, Ylower, lw=.3, color=clrs[0])\n",
    "ax.fill_between(x, Ylower, Yupper, color=clrs[1], alpha=.35)\n",
    "\n",
    "Ymean = [np.mean(yy[j]) for j in range(len(x)-1)]\n",
    "Ymax_x = x[1:][np.argmax(Ymean)]\n",
    "Ymax_y = Ymean[np.argmax(Ymean)]\n",
    "print(colored(\"Maximum (x): %.2f\"%x[1:][np.argmax(Ymean)],'red'))\n",
    "ax.scatter([Ymax_x],[Ymax_y],color='red',marker='x',zorder=5)\n",
    "Ymean_value = np.sum([x*y for x,y in zip(x[1:],Ymean)])\n",
    "print(colored(\"Mean (--): %.2f\"%Ymean_value,'green'))\n",
    "ax.axvline(Ymean_value, color='C2', ls='dashed')\n",
    "\n",
    "ax.set_xlabel('Serial interval (days)'); ax.set_ylabel('Probability density')\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "ax.set_ylim(bottom=0, top=.003)\n",
    "ax.set_xlim(0, 20)\n",
    "ytks = [0,.001,.002,.003]; ax.set_yticks(ytks); \n",
    "ax.spines['left'].set_bounds(0,.003)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Median can be calculated as $e^\\mu$, where $\\mu$ is the meanlog (param1 in our notation) of the lognormal distribution ([wikipedia](https://en.wikipedia.org/wiki/Log-normal_distribution)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34mMedian: 3.96\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "print(colored(\"Median: %.2f\"%np.exp(np.mean(param1)),'blue'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which coincides with the value stated in our paper."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}