notebooks/exponential_contact_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import integrate\n",
    "import os\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'C:\\\\work\\\\oxford_diss'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Contact Rate Model\n",
    "\n",
    "Comparing simulation results of CTMC SIR model with contact rate\n",
    "$$\n",
    "B(t) = A + Be^{-\\alpha t}.\n",
    "$$\n",
    "We compare the simulations to our novel ODE approximation model given by the solution to \n",
    "$$\n",
    "\\frac{dq}{dt} = -q^2 + (B(t) + \\mu)q - \\mu.\n",
    "$$\n",
    "In this application, we consider the contact rate function $B(t)$ with parameter values $A = 1.4, B = 1, \\alpha = 0.5$ and removal parameter $\\mu = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters of the model\n",
    "\n",
    "A = 1.4\n",
    "B = 1\n",
    "alpha = 1\n",
    "mu = 1\n",
    "\n",
    "def rate(t, A, B, alpha):\n",
    "    # contact rate for the model\n",
    "    return A + B*np.exp(-alpha*t)\n",
    "\n",
    "def exponential_decay_rate(t, A, a):\n",
    "    return A*np.exp(-a*t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating the CTMC Exponential Model\n",
    "\n",
    "We note that for a model with an exponential propensity function, instead of solving the integral for the next reaction time, we can instead solve the Lambert Function. This is described further in the paper [1]. We extend the work of the paper given above by considering a propensity function with both time dependent and time dependent parts. If we let $i$ denote an infection occuring and $r$ a removal with propensities\n",
    "\\begin{align}\n",
    "    &a_i(t) = A + Be^{-\\alpha t} \\\\\n",
    "    &a_r = \\mu.\n",
    "\\end{align}\n",
    "The probability that the reaction $i$ occurs at time $t + \\tau$ is then given by\n",
    "\\begin{align}\n",
    "    \\mathbb{P}(i, \\tau) &= (A + Be^{-\\alpha ( t+ \\tau)})\\exp\\left\\{ -\\int_0^\\tau a_i(t + \\tau') + a_r d\\tau'  \\right\\} \\\\\n",
    "      &= (A + Be^{-\\alpha ( t+ \\tau)})\\exp\\left\\{ - \\tau (A + \\mu) - f_{\\alpha}(\\tau)Be^{-\\alpha t}\\right\\}. \n",
    "\\end{align}\n",
    "where $f_{\\alpha}(\\tau)$ is the function described in by Equation (9) of [1]. That is\n",
    "\\begin{equation}\n",
    "    f_{\\alpha}(\\tau) = \\frac{1}{\\alpha}(1 - e^{-\\alpha \\tau}).\n",
    "\\end{equation}\n",
    "\n",
    "From [1], the next reaction time is therefore a solution of the following\n",
    "\\begin{equation}\n",
    "u_1 = \\exp\\left\\{  -f_{\\alpha}(\\tau)Be^{-\\alpha t} - \\tau(A + \\mu) \\right\\}\n",
    "\\end{equation}\n",
    "where $u_1$ is a randomly generated uniform number. The time to the next reaction $\\tau$ is therefore given by a solution to the Lambert function\n",
    "\n",
    "\\begin{equation}\n",
    "    \\tau = \\frac{1}{\\alpha} W\\left(  D\\exp[ D + \\beta]   \\right) - D - \\beta\n",
    "\\end{equation}\n",
    "where $D = \\frac{Be^{-\\alpha t}}{A + \\mu}$ and $\\beta = \\alpha \\log(u_1) / (A + \\mu)$. In this part, we have omitted the population weights of the propensities which must be taken into account when implementing (i.e $\\mu$ is actually $\\mu I(t)$ and $A + Be^{-\\alpha t}$ represents $(A + Be^{-\\alpha t})S(t)I(t)$ ).\n",
    "\n",
    "\n",
    "[1] - https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1334982"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def time_dependent_gillespie_algo(t0, t_final, X_0, parameters):\n",
    "    from scipy.special import lambertw\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    Temporal Gillespie Algorithm for a typical SIR CTMC model:\n",
    "    S + I -> I + I (with rate B(t) = A + B*e^{-a*t})\n",
    "    I -> R (with constant rate gamma)\n",
    "    \n",
    "    Parameters: \n",
    "    - t0: initial infective time\n",
    "    - t_final: final simulation time\n",
    "    - X_0: initial breakdown of population: (S, I, R)(0)\n",
    "    - parameters: [A, B, alpha, mu]\n",
    "    \n",
    "    Returns:\n",
    "    - t_hist: history of when reactions occur\n",
    "    - state_hist: the evolution of the epidemic.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    def calc_a_i(X, t, parameters, N): \n",
    "        A = parameters[0]*X[0]*X[1] # A\n",
    "        B = parameters[1]*X[0]*X[1]*np.exp(-parameters[2]*t) # Be^{-\\alpha t}\n",
    "        return A/N, B/N\n",
    "        \n",
    "    def calc_a_r(X, t, parameters):\n",
    "        return parameters[3]*X[1] # \\mu\n",
    "    \n",
    "    def f(tau, parameters):\n",
    "        return (1 / parameters[3])*(1 - np.exp(-parameters[3] * tau))\n",
    "    \n",
    "    def calc_a0(X, t, parameters, N):\n",
    "        A, B = calc_a_i(X, t, parameters, N)\n",
    "        A = A*X[0]*X[1]\n",
    "        B = B*X[0]*X[1]\n",
    "        mu = calc_a_r(X, t, parameters)\n",
    "        mu = mu*X[1]\n",
    "        a_list = [(A + B)/N, mu]\n",
    "        a0 = np.sum(a_list)\n",
    "        return a0, a_list\n",
    "\n",
    "    def v1():\n",
    "        # transmition, S + I -> I + I\n",
    "        return np.array((-1, +1, 0))\n",
    "    \n",
    "    def v2():\n",
    "        # recovery, I -> R\n",
    "        return np.array((0, -1, +1))\n",
    "\n",
    "    \n",
    "    S = X_0[0]\n",
    "    I = X_0[1]\n",
    "    R = X_0[2]\n",
    "    N = np.sum(X_0) # population size\n",
    "    \n",
    "    t = t0\n",
    "    state_hist = []\n",
    "    t_hist = []\n",
    "    \n",
    "    stch_vecs = {0:v1, 1:v2}\n",
    "    \n",
    "    X = np.asarray(X_0)\n",
    "    \n",
    "    while t <= t_final:\n",
    "        # calculate a_0\n",
    "        \n",
    "        A, B = calc_a_i(X, t, parameters, N)\n",
    "        mu = calc_a_r(X, t, parameters)\n",
    "        \n",
    "        # generate a reaction time\n",
    "        r1 = np.random.rand()\n",
    "        \n",
    "\n",
    "        # generate time using Lambert function (see: https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1334982)\n",
    "        alpha = B / (A + mu)\n",
    "        beta = parameters[3]*np.log(r1) / (A + mu)\n",
    "        tau = (1 / parameters[3])*(np.real(lambertw(alpha*math.exp(alpha+ beta))) -alpha - beta)\n",
    "        t += tau\n",
    "        r2 = np.random.rand()\n",
    "\n",
    "        # choose which reaction occurs\n",
    "        \n",
    "        a0, a_list = calc_a0(X, t, parameters, N)\n",
    "        for idx in range(len(a_list)):\n",
    "            if r2*a0 < np.sum(a_list[:idx +1]):\n",
    "                reaction_idx = idx\n",
    "                break\n",
    "            elif idx == (len(a_list) - 1):\n",
    "                reaction_idx = idx\n",
    "\n",
    "        # update state vector\n",
    "        X = X + stch_vecs[reaction_idx]()\n",
    "        state_hist.append(X)\n",
    "        t_hist.append(t)\n",
    "            \n",
    "    return t_hist, state_hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:78: RuntimeWarning: invalid value encountered in double_scalars\n",
      "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8576152324676514\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xV9f348df7jtybvQkBwpKwlxhBwS0g4F5VWy0qLfartVqtilqrddJqa7Wiwk8tSuugjroHKk4EDUM2MmQEAoEEsscdn98f50QCBJJAkpPxfj4e93HO+ZzPOed9IPd9PvdzlhhjUEop1T64nA5AKaVU89Gkr5RS7YgmfaWUakc06SulVDuiSV8ppdoRj9MBHEpKSorp3r2702EopVSrsnDhwl3GmNTa5rXopN+9e3eys7OdDkMppVoVEdl0sHnavaOUUu2IJn2llGpHNOkrpVQ70qL79JVSygmBQICcnBwqKiqcDuWQ/H4/Xbp0wev11nsZTfpKKbWfnJwcYmNj6d69OyLidDi1MsaQn59PTk4OPXr0qPdy9ereEZGNIrJMRJaISLZdliQic0RkrT1MtMtFRB4XkXUislREhtVYz0S7/loRmdjAfVRKqWZRUVFBcnJyi034ACJCcnJyg3+NNKRP/1RjzFBjTJY9PQX4xBiTCXxiTwOMBzLtz2TgKTvAJOBuYAQwHLi7+kChlFItTUtO+NUOJ8Yj6d45FzjFHn8e+Ay4zS5/wVjPbJ4vIgkikm7XnWOMKbCDnQOMA146ghhqtaOogn/P34TP4yLC48LncdtDF7F+70/lXrcLj0twuwSPW3CL4HG7fppfvUyE29Uq/gCUUqou9U36BvhIRAww3RgzA0gzxuQCGGNyRaSDXbczsKXGsjl22cHK9yEik7F+IdC1a9cG7MpeuYUVPDF3HY35qoAIjwuf24XXPghYBw0hwj6gRLjFHlrzon0e4vxeYv0eYv0eItwuPG5rGbfLRUKklw5xPtLi/HSI9elBRSl1gAceeIAXX3wRt9uNy+Vi+vTpjBgx4ojWWd+kP8oYs81O7HNEZPUh6taWvcwhyvctsA4oMwCysrIOK20PzUhgw4MTCIQMVaEwVcEwlcEQlYEwRRUBqoJWWSBsCIXDhMIQCocJhg3BkNlbPxiu8QlZy9jrs8YNlcEwVaEwgWCYikCYovIglcEQpZUhiioClFQG6zz4pMX5SIyKID7SS1J0BKmxPromRdE9OZreabFkJEXqQUGpduabb77hnXfeYdGiRfh8Pnbt2kVVVdURr7deSd8Ys80e5onIG1h98jtEJN1u5acDeXb1HCCjxuJdgG12+Sn7lX92RNEfgogQ4bFa3/iaait1C4cNpVVBAiFDMGQdaIKhMLvLAuQVVbBldzkrthVSUhFkT3mAdXklfLVuF8UVwZ/WkRDlZXCXBPqkxZCZFsuQLglkdojB5dIDgVJtVW5uLikpKfh8VgJLSUlplPXWmfRFJBpwGWOK7fGxwL3AW8BEYKo9fNNe5C3gtyLyMtZJ20L7wPAh8GCNk7djgdsbZS9aMJdLiPUfeA1tt+SDL2OMYU9ZgHU7S/hhRzHLcgr5PqeQBRvyqQyGAYjze8jqnsSx3ZM4uXcqfTvG6kFAqSbw57dXsHJbUaOus3+nOO4+e8Ah64wdO5Z7772X3r17M3r0aC655BJOPvnkI952fVr6acAbdveCB3jRGPOBiHwHzBaRScBm4GK7/nvABGAdUAZcBWCMKRCR+4Dv7Hr3Vp/UVfsSERKjIzg22krq2F14obBhY34pizfvIXtjAd9tLODT1Xn85YPV9O0Yy8VZGZyYmUJmhxjtDlKqlYuJiWHhwoV8+eWXzJ07l0suuYSpU6dy5ZVXHtF6pSW/GD0rK8voUzYPbeuecj5asZ3/LNjMurwSAEb368Dlx3XjuJ7J+L1uhyNUqvVZtWoV/fr1czqMfbz66qs8//zzvP322/uU1xariCyscXn9PvSO3Fauc0IkV43qwZUju7O5oIzXF23lqc/X8/GqPLolR3HnhH6M6Z+mLX+lWpk1a9bgcrnIzMwEYMmSJXTr1u2I16tJv40QEbolR/P7Mb35zclH8fkPeTz84Romz1pI346x/H5Mb07unaotf6VaiZKSEq6//nr27NmDx+OhV69ezJgx44jXq0m/DYqMcDNuYDqn90vjxQWbeWLuOq6ZtZA4v4cLj+nCdaf2IiXGwUualFJ1OuaYY5g3b16jr1eTfhvmdbuYOLI7P8vKIHtTAbOzc3h+3kZe/nYLN4/tzVmDO9Ex3u90mEqpZqTP028HIiPcnJiZyj8vO5r3bziJgZ3juP/dVZzyyFz+s2AToXDLPZmvlGpcmvTbmT4dY5l9zfG8fu1I+naM4843lnPBU/PYWVzpdGhKqWagSb8dEhGGdU3ktf8bydQLBrFmexHjH/uCeet2OR2aUqqJadJvx9wu4dLhXXnj2lEkRkXw82cW8NsXF1FcEXA6NKVUE9Gkr+iXHsdr147k0mMzeGdpLudN+5oNO0ucDksp1QQ06SsA4vxepl44mOeuzCK/tIpLZsxn7Y5ip8NSqt2KiYlpkvVq0lf7OK1vGv+95niMgZ8/s4DcwnKnQ1JKNSJN+uoAmWmxPH35MIorAlz01DdsL2zYOziVUi2X3pylapXVPYnnrxrO1TO/46Kn5/Hir46ja3KU02Ep1fzenwLblzXuOjsOgvFTG3ed9aQtfXVQI3om8+Kvj6O4IsjkWdlU2c/yV0q1XtrSV4c0JCOBv140mGtmLeSWV7/nkYuH4HVrW0G1Iw61yJuKfntVnc4Y0JGbx/TmzSXbePqz9U6Ho5Q6AtrSV/Xy29N6sXpHMX+b8wPBsOH3Y3o7HZJS6jBoS1/Vi4jw958N4YReKTwxdx1Lc/Y4HZJSbVpJSdPcIKlJX9Wbz+PmH5cOJS3Wx9Uzs9lSUOZ0SEqpBtKkrxokJcbHzKuHU1IZ4JZXvycQ0it6lGpNNOmrBuudFsu95w5k/oYCHnh3ldPhKKUaQJO+Oiw/y8rgl8d3Y+a8jXz+w06nw1FK1ZMmfXXY7pjQj54p0dzx+jK9cUupVkKTvjpsfq+bu88ZwNY95cyav8npcJRS9aBJXx2Rk3uncmz3RGZ8sV5fvqJUI3K73QwdOpSBAwdy9tlns2dP41wmrUlfHbE7JvQjr7iS//flj06HolSbERkZyZIlS1i+fDlJSUlMmzatUdarSV8dsaO7JjK6Xxr//HQt7y7NdTocpdqc448/nq1btzbKuvQxDKpRPHTBIPKer+T3s5fQLTmKgZ3jnQ5JqUbxl2//wuqC1Y26zr5Jfblt+G31qhsKhfjkk0+YNGlSo2xbW/qqUaTE+HhuYhYp0RFcM2shhWXav6/UkSgvL2fo0KEkJydTUFDAmDFjGmW92tJXjSY5xseTlx/D+U9+zT8++YG7zx7gdEhKHbH6tsgbW3WffmFhIWeddRbTpk3jd7/73RGvV1v6qlENzUjgkqwMZn2ziTXb9cXqSh2p+Ph4Hn/8cR555BECgSP/BV3vpC8ibhFZLCLv2NM9RGSBiKwVkVdEJMIu99nT6+z53Wus43a7fI2InHHE0asW6eaxffC4hQffW4UxxulwlGr1jj76aIYMGcLLL798xOtqSEv/BqDmg1b+AjxqjMkEdgPVZxkmAbuNMb2AR+16iEh/4FJgADAOeFJE3EcWvmqJUmN9/GFsHz7/YSfvLtOreZQ6HPs/Wvntt9/miiuuOOL11ivpi0gX4EzgGXtagNOAV+0qzwPn2ePn2tPY80+3658LvGyMqTTG/AisA4Yf8R6oFumqUT0Y3CWee95aye7SKqfDUUrZ6tvS/wdwK1D9gJVkYI8xJmhP5wCd7fHOwBYAe36hXf+n8lqW+YmITBaRbBHJ3rlTH+TVWrldwkMXDCK/tJL/9+UGp8NRStnqTPoichaQZ4xZWLO4lqqmjnmHWmZvgTEzjDFZxpis1NTUusJTLdiATvGcNbgTz371o75wRbU6reF81OHEWJ+W/ijgHBHZCLyM1a3zDyBBRKov+ewCbLPHc4AMAHt+PFBQs7yWZVQbdceEvnhcws2zvycUbvlfIqUA/H4/+fn5LTrxG2PIz8/H7/c3aLk6r9M3xtwO3A4gIqcAfzDG/EJE/gtchHUgmAi8aS/ylj39jT3/U2OMEZG3gBdF5O9AJyAT+LZB0apWJz0+krvPHsCtry3l7e+3cd7RB/ToKdXidOnShZycHFp6F7Pf76dLly4NWuZIbs66DXhZRO4HFgPP2uXPArNEZB1WC/9SAGPMChGZDawEgsB1xpjQEWxftRIXHdOF6V+s56nP1nPm4HS8br09RLVsXq+XHj16OB1Gk5CW/PMlKyvLZGdnOx2GagTvLN3Gb19czMTju/Hncwc6HY5SbZqILDTGZNU2T5tcqlmcOSidy4Zn8ML8TSzfWuh0OEq1W5r0VbMQEaaM70dSVAR3v7WCsJ7UVcoRmvRVs4mP9DJlfF8WbtrNC99sdDocpdolTfqqWV10TBdO6JXC3z76Qa/dV8oBmvRVsxIRHjx/EAhMnrWQsqpg3QsppRqNJn3V7LomR/GPS4ayKreI2d9tqXsBpVSj0aSvHHF6vzSGZiTw9OcbKKrQt2wp1Vw06SvH/PmcAeQVV/DAO6vqrqyUahSa9JVjhmQkMPmko3glewvz1u1yOhyl2gVN+spRN47OJCUmgr98sJrKoD6VQ6mmpklfOcrvdfOnswfwfU4h0z5d53Q4SrV5mvSV484Z0okLju7ME3PX8eXalv1UQ6VaO036qkW488x+pMdHcu/bK1v0M8yVau006asWITnGx42jM1mbV8Ibi7c6HY5SbZYmfdVinHd0Z4ZkJPDAu6soLNNr95VqCpr0VYvhdbt48PyB7C6r4q8frnY6HKXaJE36qkUZ0CmeK0f24MVvN/PtjwVOh6NUm6NJX7U4N43tTbekKH730mJ2l1Y5HY5SbYomfdXixPg8PPHzYRSUVnHT7CX6whWlGpEmfdUiDewcz23j+zJ3zU5eW5TjdDhKtRma9FWLddXI7vTtGMvDH65hV0ml0+Eo1SZo0lctlsslTL1wMHvKAtz26lJC2s2j1BHTpK9atKEZCdw6rg+frM7jybn6bB6ljpQmfdXiXT2qB6f2SeXRj3/gwxXbnQ5HqVZNk75q8Vwu4bHLjqZnagzXv7iYNduLnQ5JqVZLk75qFeL8Xl6ZfBwxfg+3vPo9wVDY6ZCUapU06atWIznGx33nDmRpTiHTv9jgdDhKtUqa9FWrcubgdM4clM5jH69l465Sp8NRqtXRpK9anbvO6o/P4+IP/9VuHqUaSpO+anU6xvu548x+ZG/azaz5m5wOR6lWpc6kLyJ+EflWRL4XkRUi8me7vIeILBCRtSLyiohE2OU+e3qdPb97jXXdbpevEZEzmmqnVNt3SVYGJ/VOZer7q1m/s8TpcJRqNerT0q8ETjPGDAGGAuNE5DjgL8CjxphMYDcwya4/CdhtjOkFPGrXQ0T6A5cCA4BxwJMi4m7MnVHth8slPHLxYCI8Lm54eTFlVUGnQ1KqVagz6RtLdVPKa38McBrwql3+PHCePX6uPY09/3QREbv8ZWNMpTHmR2AdMLxR9kK1Sx1i/Tx80RBWbCvi7jdX6Lt1laoHT30q2S3yhUAvYBqwHthjjKluXuUAne3xzsAWAGNMUEQKgWS7fH6N1dZcpua2JgOTAbp27drA3VHtzbiBHZl4fHdmztuIxy08eP4grDaGaouMMVSGKqkKV1EVqiIQCvw0XhW2p+3xmmXBcJCwCRMmjDEGYwxhwoSNNR02YQwHlhvM3joHK7eXrTkMG+sCg+rxmvVr1qleb23lw9KGcdXAqxr937BeSd8YEwKGikgC8AbQr7Zq9rC2b5w5RPn+25oBzADIysrSppuq05/O6g/AzHkb6ZkSw69P6ulwRGp/xhiKA8UUlBdQHiynIlRBRdD6VIYqKQ+Wk1+Rz/bS7RRUFBAIW4k6GA5SEaygsKqQwspCiiqLCBrnuvIEwSUuRAQXrp/G9ykX197pGuXV44cst9crIhRWFjbJPtQr6VczxuwRkc+A44AEEfHYrf0uwDa7Wg6QAeSIiAeIBwpqlFeruYxSh83lEu4+uz/r8kp45KM19E2P5cTMVKfDahcqQ5XkFOewo2wH+eX5FFUVURWqIr88nx1lO9hVvotd5bvYXrqdilBFnetL8CWQ7E/G6/biEQ8elwefx0dmVCbxvngSfAlEeaOIcEUQ4bY+XpfXGj9Emcfl2SfhVifn6gS7fxLfp16N8ragzqQvIqlAwE74kcBorJOzc4GLgJeBicCb9iJv2dPf2PM/NcYYEXkLeFFE/g50AjKBbxt5f1Q7JSL8/ZIhXPjUPH730mL+d90ouiVHOx1WmxAIBdhetp0fC3/kx8If2Vi0kc1Fm9lcvJkdpTswB/5gx+vy0immE8n+ZPok9eGEzifQMbojSf4kor3R+N1+/B4/Po+PSHckPo+PRF8iUd4oB/awfZG6Tn6JyGCsE7NurBO/s40x94pIT6yEnwQsBi43xlSKiB+YBRyN1cK/1BizwV7XncDVQBC40Rjz/qG2nZWVZbKzs49k/1Q7s3ZHMRc8OY/4KC//uvJYMtNinQ6p1SgNlLJ813Kyd2SzrWQbuaW5bC/dTm5pLsHw3i6VBF8C3eK60S2uGxmxGWTEZvyU0BN9iUS4I/C5fbhdenGeU0RkoTEmq9Z5LfmKB0366nAs3LSbyS9kkxrr4/VrRxIV0aBezDavNFDKhj0byCvLY92edawvXM+q/FVsLNoIgEtcdIjqQHp0Oh2jOpIek073uO50j+9Ot7huJPmTnN0BVadDJX39Nqg255huiTx88WB+9Xw2v39lCU/+4hjcrrbRH3s4SqpKWLB9AasLVjNv6zxW5K8gZEI/zU+LSqN/cn8m9JxAn8Q+DO84nJiIGAcjVk1Jk75qk07rm8Zt4/ry0PurueP1ZUy9sH1dyllQUcDXW7/m3R/fJXt7NpUh6x3DA5MHcvXAqxmUMohEfyJ9k/ri9/gdjlY1J036qs2afFJPCkqrmP7FBrqlRHHtKb2cDqnJhMIhFuQu4MutXzJ3y1y2lmwFoENkBy7qfRFjuo0hMzGTuIg4hyNVTtOkr9osEeG2cX1Zv7OEv36whozEKM4e0snpsBrVxsKNzP5hNh9v+pjc0lw8Lg8ndDqBS/tcysCUgQxOHUyEO8LpMFULoklftWkul/C3nw3lyn99y82zv6dHSjQDO8c7HdYRKQ2U8v3O73l59ct8nvM5xhhGdR7FTcfcxCkZp2h3jTokvXpHtQu5heVc8OQ8isoDPHrJUMYO6Oh0SA1SEaxg7pa5vLvhXb7e9jXBcJBITyS/6PcLLu1zKWnRaU6HqFoQvWRTKazE/8tnvyVndzkzrzqWET2TnQ6pTmsK1vDKmlf44McPKA4U0yGqA+O7j+eYtGMYkT5Cb2ZStdKkr5Qtr7iCnz39DduLKrj3nIGcPaQTkREt6yaiQDjA3M1z+e8P/2V+7nz8bj+ndzud83udT1Zalt70pOqkSV+pGvKKKvjFMwtYm1fCkIwEnvzFMDonRDodFkVVRby57k2eW/4cu8p3kR6dzkW9L+LSvpfqVTeqQTTpK7WfcNjwSvYW7vrfcvxeN69fO5LeDj2yobCykH+v+jevrH6F3ZW7GZwymGuGXMOoTqO0Va8Oi96Rq9R+XC7hsuFdGdQ5nqtmfscFT85j2i+GcXLv5n0655aiLVz7ybVsLNrIyE4juW7odQxOHdysMaj2RV+Mrtq1gZ3jefO6UWQkRXH1zO945bvNzbbt7O3ZXPLOJeRX5PPs2GeZPma6JnzV5DTpq3avU0Ik//3N8YzqlcJtry3jjjeWUVwRaLLtVYWqeG75c0z6aBJJkUm8OOFFhqfrm0NV89CkrxQQ4/Pw7MQsJp3Qg5e+3czov3/O8q2N/+airSVbuezdy3h04aMM7zicl898me7x3Rt9O0odjCZ9pWxet4u7zurPK5OPpzIY5ppZC1m7o7hR1l0ZquS1H17jwrcuZEvxFu4YcQczxszQp1mqZqdJX6n9DO+RxHNXHktRRYCz/vkV0+auozIYqnvBg/jgxw8449UzuOebezgq4SheO+c1Lut7Wbt66qdqOTTpK1WLYV0T+ej3JzG8RxIPf7iGc/75NevyShq0jkA4wEMLHuKWL26hY3RHZoyZwQvjXiAjNqPuhZVqIpr0lTqI9PhIXrh6OM9OzGJbYTkXPT2PZ77cQCAUrnPZQDjAlC+m8OLqF7m83+XMmjCL4zsdr9fdK8dp0lfqEESE0/ul8ca1I0mN8XH/u6u48Kl5bC+sOOgy20q28Zs5v+GjTR9x47AbuW34bXhd3maMWqmD0ztylaqncNjwvyVbmfL6MtwiZHVP5HenZ5LVLfGn/vkV+Su47uPrKA+WM2X4FM7PPN/hqFV7pHfkKtUIXC7hgmFdyOwQy0vfbebdpblc/PQ3TBjUkYcuGMyCHXP507w/ER8Rz3NnPEfPhJ5Oh6zUAbSlr9RhKiwL8K95P/L4pyuJ7fIm4ehsMqIzmT52Ghlx6U6Hp9oxbekr1QTio7xMPCGVL4pns7ZoGZW7TmHl6tM5b9VyjuuZy8VZXTgpMxWXSy/NVC2HJn2lDlNuSS6TPprEjtIdPHzSwxybeirz1ufzv8Vb+WxNHu8uyyUlJoKRR6UwuEs8ZwzoSEaSvvREOUu7d5Q6DGt3r+X6T6+nsLKQp0Y/xdAOQ/eZXxkM8e7SXD5asYNlWwvZuqccETi+ZzJnDk5n/MB0kqL1heWqaejz9JVqRBsKNzDx/Yl4XV4eOfkRhqUNq3OZjbtK+d+Srcz6ZhP5pVX4PC6O7Z5EVvdEzhyUTq8OMXqHrmo0mvSVaiS5Jblc9u5lhEyI58c/T8/4hl2hY4xh9fZiXvluC1+u3cmGXaUYA0nREfROiyE5xkesz0P3lGhO6JVC//Q4PSegGkxP5CrVCKpCVUz5cgplwTJenPBigxM+WDd79UuP455zBgCwvbCCOSu3s2xrIRt2lrJqWxFFFUF2lVQCkBwdwfFHJdOrQwx9O8aREhNB58RI0uOdf72jap006StVD2ET5v7597MobxEPnfgQvRJ7Ncp6O8b7ueL47geU7yiq4Ku1u/hq3S6+WZ/PO0tzf5onAj1SoumZEk3/9DjSEyLp0zGWo1JiiIv0aDeROiRN+krVw4ylM3hj3Rtc1Psizup5VpNvLy3Oz4XHdOHCY7oAUBEIsTK3iOKKIAs37WZ1bhGLNu/h41V5+ywX6/fQq0MMHeP89EyNZlDnBLokRuL3ukmP9xPt0698e6d/AUodgjGG6UunM23JNM7ofgZ3HXeXI3H4vW6GdU0E2Oc9voFQmG17ypm/IZ/iiiAb80tZn1fKmh3FfLRyB6Hwvufs0uJ8REV4iIpwkxQdQXJ0BMkxPpKiI0iJiSAxKsKa73MTFeHG63bh87iIi/QS69NfEW1BnUlfRDKAF4COQBiYYYx5TESSgFeA7sBG4GfGmN1i/VU8BkwAyoArjTGL7HVNBP5or/p+Y8zzjbs7SjWup79/mie/f5IJPSbw4AkP4pKW9YxCr9tFt+RouiVHHzCvIhBiVW4Ru0qqKKsKsmFnKbmF5VQEwpRUBikorWJjfikFJVWUVtX9vgCvW+wDhY+0OB+psT4iPC78Hjcd4/3E+Dy4XYLX7cLtEjwuITLCTVSEB7/XhUvE+riwx61zHNXjLhF8Hhcxfg8R9jr0INP46tPSDwI3G2MWiUgssFBE5gBXAp8YY6aKyBRgCnAbMB7ItD8jgKeAEfZB4m4gCzD2et4yxuxu7J1SqjH8LftvzFwxkzN7nskDox5odY9F9nvdHG3/OqhLRSBEfmkVu0urKA+EKKkMUl4VoioYpjIYoqg8SH5pFQWlleSXVLGpoIzV24upCoYpqwpRHjj8l8wcjAhEuF1EuF14PdVD66BSfVCoWfencQ4s3+fQUaOyHFi0T93qg86+ZbVvi1q2VVvdfcoOMf+4nslcd2rjnDuqqc6kb4zJBXLt8WIRWQV0Bs4FTrGrPQ98hpX0zwVeMNa1oPNFJEFE0u26c4wxBQD2gWMc8FIj7o9SjeKdDe8wc8VMxnYb2yoTfkP5vW46J0TSOaHhVwUZYygsD1BWFSIUNgTDhlA4TFXQUB4IUV4VoiIQImwMYWPVDxvsaYP5aRzrgFMRJBAKEwiFqQqFqQpa44GgsabtsurLzWtedW72i+vAsgPrHuyy9epiU2MN+yy/z7pqj+WnaVNd6yDrqiWWiiY4kEID+/RFpDtwNLAASLMPCBhjckWkg12tM7ClxmI5dtnByvffxmRgMkDXrl0bEp5SjeLTzZ/yx6/+yODUwTx44oNtPuEfKREhISqCBH3CRKtQ7w5KEYkBXgNuNMYUHapqLWXmEOX7FhgzwxiTZYzJSk1NrWURpZrOexve48a5N5KZmMn00dPxuX1Oh6RUo6pX0hcRL1bC/48x5nW7eIfdbYM9rL52LAeo+RLQLsC2Q5Qr1SIs3bmUO7++k6EdhjJz3ExiImKcDkmpRldn0revxnkWWGWM+XuNWW8BE+3xicCbNcp/KZbjgEK7G+hDYKyIJIpIIjDWLlPKcUVVRdz6xa2kRKbwz9P+SbT3wKthlGoL6tOnPwq4AlgmIkvssjuAqcBsEZkEbAYutue9h3W55jqsSzavAjDGFIjIfcB3dr17q0/qKuUkYwx3fnUn20u3869x/yLeF+90SEo1mfpcvfMVtffHA5xeS30DXHeQdT0HPNeQAJVqas8se4bPtnzGDcNu4OgORzsdjlJNqmXdaaJUM/si5wseX/w443uMZ9LASU6Ho1ST06Sv2q2V+Su5Ye4NHBV/FPccf4/e/anaBU36ql0KhALcP/9+Yr2xzBw3kyivXmSu2gdN+qrdCZswf/z6jyzbtYwpw6eQ4E9wOiSlmo0mfdXu/Hvlv3nvx/e4YdgNTOg5welwlGpWmvRVu7KmYA2PLXqMk7ucrCduVbukSV+1G8YYHlzwIH6Pn3tG6olb1T5p0lftxsZJixYAABObSURBVIebPmRR3iJuPOZGUiJTnA5HKUdo0lftwvbS7Tw4/0H6JPbhgl4XOB2OUo7RpK/avPzyfCa+P5GyYJk+Klm1e/qOXNWmBcIBbvniFvLK85h22jR6J/Z2OiSlHKUtfdVmGWP409d/4rvt3/HHEX9kZOeRToeklOM06as2a9bKWbyz4R2uGnAVF/a+0OlwlGoRNOmrNmll/koeW/QYozqP4sZjbnQ6HKVaDE36qs35YfcPXPvxtSRFJvHQCQ/hEv0zV6qafhtUmxIKh7hx7o24xMX00dNJ9Cc6HZJSLYomfdWmzFg6gy3FW7jl2FvomdDT6XCUanE06as24+NNH/Pk908yqtMoxnYb63Q4SrVImvRVm7AkbwlTvpzCoJRBPHrqo3oDllIHoUlftXrBcJA/f/NnkvxJPHH6E0R6Ip0OSakWS+/IVa3e1G+nsm7POv5xyj9I8ic5HY5SLZq29FWr9sbaN3hlzStc3PtiTut6mtPhKNXiadJXrdaq/FXcP/9+RqSP4I4Rd+jz8ZWqB036qlXaUbqDX8/5NTERMdw38j48Lu2pVKo+9JuiWp2wCXPrF7dSEaxg9tmzSY9JdzokpVoNbemrVmf60uksylvEnSPupGe83oClVENo0letyqr8VTz9/dOM7jqac3ud63Q4SrU6mvRVq/LCyhfwu/38edSf9UFqSh0G/daoVmNB7gLe2fAOF/e+mLiIOKfDUapV0qSvWoWKYAV//e6vdIruxPXDrnc6HKVarTqTvog8JyJ5IrK8RlmSiMwRkbX2MNEuFxF5XETWichSERlWY5mJdv21IjKxaXZHtUUFFQVM+mgSP+z+gVuOvQWf2+d0SEq1WvVp6c8Exu1XNgX4xBiTCXxiTwOMBzLtz2TgKbAOEsDdwAhgOHB39YFCqUMxxnDzZzezctdKpp44ldHdRjsdklKtWp1J3xjzBVCwX/G5wPP2+PPAeTXKXzCW+UCCiKQDZwBzjDEFxpjdwBwOPJAodYA31r1B9o5sbht+G2f2PNPpcJRq9Q63Tz/NGJMLYA872OWdgS016uXYZQcrP4CITBaRbBHJ3rlz52GGp9qCXeW7ePi7h8lKy+JnfX7mdDhKtQmNfSK3toefmEOUH1hozAxjTJYxJis1NbVRg1OtRzAc5M6v7qQyVMldx9+ll2cq1UgO95u0w+62wR7m2eU5QEaNel2AbYcoV6pWTy55knnb5nFz1s16161Sjehwk/5bQPUVOBOBN2uU/9K+iuc4oNDu/vkQGCsiifYJ3LF2mVIHWLhjIc8uf5Yze57Jz/v+3OlwlGpT6nzgmoi8BJwCpIhIDtZVOFOB2SIyCdgMXGxXfw+YAKwDyoCrAIwxBSJyH/CdXe9eY8z+J4eVoriqmJs+u4m0qDTuOu4ufVyyUo2szqRvjLnsILNOr6WuAa47yHqeA55rUHSq3fnn4n9SUFHArPGziPZGOx2OUm2Onh1TLca8rfN4afVLXN7vcoZ2GOp0OEq1SZr0VYuwp2IPd319Fz3je3LDsBucDkepNktfoqJahHvn30t+RT6Pn/Y4fo/f6XCUarO0pa8ctzhvMXM2zWHy4MkMSBngdDhKtWma9JWjlu9azjVzrqFDVAeuHHCl0+Eo1eZp0leOCYVD3Df/PrwuL9NHTyfKG3XoBUytN3ErpRpA+/SVY55b/hwr81fylxOn0ivhKAiHYM6fYOWbULjl4Av+6lP4fCp4fLDq7b3laYNgxzIYfg2c8SCU7oQ4fWm6UjWJacGtp6ysLJOdne10GKoJfJHzBTd8cj0nlJbweN6uWh/O1CSOvgLOfaK5tqaUI0RkoTEmq7Z52tJXzStYydsL/sY9P7xI70AVD+7MPzDhx3SEqGToOwEqi2HYL8GEoeMgeH8KLHjKqtf/XDh3GngiYc17MPsK8CdAxZ6Db3/xLOsz4AK48FlwaQ+nal+0pa+aR2EOPDqA7/w+ft2xA/2qqpi2fSdJ4TBc9QFkjGi6BGwMiEBVGTy4X3fPZa9AH321g2pbtKWvnPXNNPjwDgpdLm5PTSYjEGT6sFuJy/pV82y/+vk9EVFwTyHkLoXpJ1plL11iDQddDOc9DW79Sqi2Tf/CVdO6Jx6AIpfwfz36kE8lL4x7gbjUQc7FlD7YSv7v3QLfzrDKlv3X+gy+BM55AjwRzsWnVBPS7h3VNFa+CbN/CUC5CNcNOpnFJZv5y0l/YWz3sQ4HV4MxsOlrmHmQVzGe808Y9DPw6l3CqvU4VPeOJn3V+J49A7bMB2CX28Xtx5zN/J2L+UPWH5g4YGIdCzsoHIY3roFls2uff9RpcOlLegBQLZ4mfdX0Ns+H587Yp6hwwl+5ePNr7CzbyU1ZN3FF/yscCu4wBCshJxs+vQ82f7PvvN8uhJRezsSlVD0cKunr9WrqyK2ds2/C98VTetNK/m/n5+ws28kzZzzTuhI+WDd+dR8FV38Ad+XDiP/bO++JY6xzFQUbnItPqcOkLX11ZEJBuC/ZGv/5bOh6PAWhCq7/7PesyF/BQyc+xPge452NsTHtWmcl/ZrOeAgGXQQxHZyJSan9aEtfNY389XsT/gm/h95nsKJ0Cz//4JesyF/B1JOmtq2ED1a3zj2FcHWNVzx/eDs8kgn/uRjWz3UuNqXqQVv6quHKCmDW+ZC7ZG/Z3XtYt2c9l717GVHeKO4deS8nZ5zsXIzNZcNnsGgW7FgOFYVQnGuV601fykF6c5ZqPMtehdcm7Z3+5VvQ82Syt2cz5cspRHmjePXsV0mNSnUuxubU8xTrAxCogHmPw9wH9t70dcX/4KhTnYlNqVpo946qn5Kd1snLmgn/TwXs7Nifmz67ias+vAqDYfqY6e0n4e/P64eTb4WbVlvX9gPMOg/evA6Kdzgbm1I27d5RddvwGbxwrjXuj7ceV9B3Auv3rOfaj68lrzyPkzqfxD0j7yHRn+hoqC1KVSl8MAUWvbC3zOWB0/4IQ3+hJ35Vk9Hr9NXhCVTAx/fsfarlhc/CoIsIhoO8sPIFnv7+aSI9kTx26mMM7TDU0VBbtJyF8N3/g+9fOnDe/30Daf2bPybVpmnSVw237hP49wXWuDcaLn8Nuh3PivwV3PbFbWwq2sRx6cdxz8h76BzT2dlYWxNjrH7/H7+AdR9bZal9rb5/feGLaiSa9FX9lRXAX3vsnT71TjjhJjaXbuPZ5c/y+trX6RDVgVuOvYVx3fXqlCOyZwvMnAB7NlvTvcbAsb+CzLH6nH91RPTqHXVohVvh5csgWAU7V/9UXH79Qt4vWMqHn17HvG3zEISLel/E9UdfT5I/ycGA24iEDLhxGeSthv9eCes/hXVzIKEb9BoNcZ0g+SjoOBgSe+iBQDUKTfrtUaDCal2ufgfmPgjhwN5ZnYay8cQbeGnPcj75eBIFFQWkR6fzmyG/4fxe59MpppODgbdRHfrCdfOhfA+s/cg68bvoeQgH99bxxUFKpvVrIL4zdDnWOij44va+L0CpetDunfagKBfyVsKuH6y+5B+/gKqSvfN7jWFZ/3H8a88yvtw2j4pQBV6Xl1MzTuXiPhczouMIRBNL8wqHYfv3sOkb8EbC+k9gx4oDn/eTnAn9z4HUftD7DPDHOROvalG0T7+9ylsNn0+FFW/sW95lOPQaTV5cGm+YIt7b8gkbCjcQFxHH+B7jGZA8gJGdRpIWneZM3OrgirfDmvfBFwu7N1ovfqnRJUenoyE+AxK6Wk8KTT7KOhgk9XQsZNX8NOm3F8bAnk1WIlj4AhRuBnFbLxDvNpKC9IEsDhSydPcq5m6Zy4+FPwKQlZbF6V1P5/zM84n2Rju8E6pBjIGda6wXwRRssB4HUZhjjUfEQmWhVS++K/hirOSfPtQ6GABERENsulXui3FuP1Sj0qTfFhljtfp2/WD97P/+JSj4kfxgKbkeD7keN7t7nUZR12NZW7adpTuXklOSA4BHPPRP6c/orqMZ1XkUvRN7O7wzqtGFQ+Byw6611gnizd9YdwUX5liNgdq4vNYL6jOOhbjOEJkI/gTrV0NENEQlWV1NqsVrUUlfRMYBjwFu4BljzNSD1W23ST9QARV7rMsnC9bDzjVUFm4m3wTJLd3O7pJtlBbnsktC7HB72Ob1sMnnZ6fbTSnhA1aXGpnKkNQhDE4dzKCUQQxMGYjfo29/arcqi61fAi6v9XdWkgf566xfC2X5sH05mFDty7q8EJUM8V2sT1QSeKOsT0SUdU9HhD3tj4fIBHD7rLKoZD3x3ExazCWbIuIGpgFjgBzgOxF5yxizskk3bIzVvxmsgMoi66oIY6zWS0QMhAOEA+UEg5UEwkECJkjAhAiKEA5WEQ5VYcIBwiF7PFRFOBQkHK6yywKEwtYwGAoQDlcRDFURDgcIBasIBUqpCgepMiEqTYgqE6YyVEllZSFVlUVUBsspw1AWDlAeDlBCmEoRKkUodLvY4dnvvykCSLZO2MW6/XSM7kSvhB6Miu5Ip+hOdI7tTOeYziT4EoiLiCPKG9Wk/7yqlfHFQvqQWmb8wRoEq6yDQflu66mhJTutE/9l+dawdKf1i2H7MuvJooEy61Mf4rYPEn7rV4Mn0hp6I61fFjFp1vfShK1fKyZsfTx++/ta4+Pxg9sL7ghr6KoxXl3u8lgft3fvePVHalwCK2JNi7vNXxrb3JdsDgfWGWM2AIjIy8C5QKMm/R/WvM0tX00hbAxhDCHAAGGBMBBGCAsEgYAIARGCDrU+vIDP68KPi2hXFD6Xl2i3jzhxExERR7+YNGJj0gm53PRK6EWXmC4kRyYT5YkiKTJJ++BV4/NEWM8FiukAqX3qt0w4DMFyCJRbzxwKlFm/IALlEKq0ysoKoLzAKqv+BMutX7aBMti9CbYssKZdbjsR28OgvQ6aqWdC3HYM1UOX9XG5rYNR9QFJBBA71prjrrrHD1jOte945hg444FG37XmTvqdgS01pnOAETUriMhkYDJA165dD2sjkTFp9IrqiEs8uFxuXOK2hi4PLneENUTwmBDeUAiv24vXFUGEOwKvuPCIGy+CxxhcLi8utwcRDy63B5fLg7i89vq8iMuL216nxxOByxWB2+3D7fHaQx+uiDh83kgi3BH43D58bt9P4y5p260K1U64XHtb4NEpVlmHfo27DWP2HlSqiq1fJKEqCAWse01C1dPBvePhkDUvHLQ+IXsYDljr+2ndYWvaVCf0GkNjaoyH7QOBx0rQxgCmxjBcj3Gz33i49vG4pnm8SXMn/dqa0/scuo0xM4AZYPXpH85GMjoP52+XfXI4iyqlWioR69xARBTQTh/f3Qiau5mZA2TUmO4CbGvmGJRSqt1q7qT/HZApIj1EJAK4FHirmWNQSql2q1m7d4wxQRH5LfAh1iWbzxljVjRnDEop1Z41+wPXjDHvAe8193aVUkrpO3KVUqpd0aSvlFLtiCZ9pZRqRzTpK6VUO9Kin7IpIjuBTYexaAqwq5HDaW66Dy2D7kPLoPvQMN2MMbXewdaik/7hEpHsgz1hrrXQfWgZdB9aBt2HxqPdO0op1Y5o0ldKqXakrSb9GU4H0Ah0H1oG3YeWQfehkbTJPn2llFK1a6stfaWUUrXQpK+UUu1Im0r6IjJORNaIyDoRmeJ0PA0lIhkiMldEVonIChG5wemYDpeIuEVksYi843Qsh0NEEkTkVRFZbf9/HO90TA0lIr+3/46Wi8hLIuJ3Oqb6EJHnRCRPRJbXKEsSkTkistYeJjoZY10Osg8P239PS0XkDRFJcCK2NpP0a7x0fTzQH7hMRPo7G1WDBYGbjTH9gOOA61rhPlS7AVjldBBH4DHgA2NMX2AIrWxfRKQz8DsgyxgzEOtR5pc6G1W9zQTG7Vc2BfjEGJMJfGJPt2QzOXAf5gADjTGDgR+A25s7KGhDSZ8aL103xlQB1S9dbzWMMbnGmEX2eDFWommaF2U2IRHpApwJPON0LIdDROKAk4BnAYwxVcaYPc5GdVg8QKSIeIAoWslb6owxXwAF+xWfCzxvjz8PnNesQTVQbftgjPnIGBO0J+djvTmw2bWlpF/bS9dbXcKsJiLdgaOBBc5Gclj+AdwKhJ0O5DD1BHYC/7K7qJ4RkWing2oIY8xW4BFgM5ALFBpjPnI2qiOSZozJBatxBHRwOJ4jdTXwvhMbbktJv86XrrcWIhIDvAbcaIwpcjqehhCRs4A8Y8xCp2M5Ah5gGPCUMeZooJSW352wD7vP+1ygB9AJiBaRy52NSgGIyJ1YXbn/cWL7bSnpt4mXrouIFyvh/8cY87rT8RyGUcA5IrIRq4vtNBH5t7MhNVgOkGOMqf6V9SrWQaA1GQ38aIzZaYwJAK8DIx2O6UjsEJF0AHuY53A8h0VEJgJnAb8wDt0k1ZaSfqt/6bqICFY/8ipjzN+djudwGGNuN8Z0McZ0x/o/+NQY06pamMaY7cAWEeljF50OrHQwpMOxGThORKLsv6vTaWUno/fzFjDRHp8IvOlgLIdFRMYBtwHnGGPKnIqjzSR9+wRJ9UvXVwGzW+FL10cBV2C1jpfYnwlOB9VOXQ/8R0SWAkOBBx2Op0HsXymvAouAZVjf9RbxGIC6iMhLwDdAHxHJEZFJwFRgjIisBcbY0y3WQfbhCSAWmGN/t592JDZ9DINSSrUfbaalr5RSqm6a9JVSqh3RpK+UUu2IJn2llGpHNOkrpVQ7oklfKaXaEU36SinVjvx/RhDhcgoS1q4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "start = time.time()\n",
    "t_hist, state_hist = time_dependent_gillespie_algo(0, 50, [4999, 1, 0], [0.7, 2, 0.02, 1])\n",
    "end = time.time()\n",
    "print(end - start)\n",
    "plt.plot(t_hist, state_hist)\n",
    "plt.legend([\"S\", \"I\", \"R\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAErCAYAAADEyxRmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXwV5dn/8c9FwhJCQkD2NQiCoihiEBQX6oK71gWXat2lVq3a5ana6s/ap62PrVqX1rpLtUqt+1ZxxxVFQGQRRECWyL6FLQRCrt8fM4mHcJJM4CxZvu/Xa17nnJn7zFxnzpxcmbnvuW9zd0RERCprku4ARESkblKCEBGRuJQgREQkLiUIERGJSwlCRETiUoIQEZG4lCBERCQuJQgREYlLCUJEROJSghARkbiUIEREJC4lCBERiUsJQkRE4lKCEAHM7Bgz+zBCuafN7JJUxCTpYWbjzOzGdMdRFyhBpIGZFZjZi2a2wszWmdlsM7vLzDonaP0JPcCjrC8sU2JmG8ysyMymmNnIRG9nZ5nZYjMrDuMrNrO5ZnZRuMyAvwI3x3nfm2b2i5hZNwN/MrOsZMRZ3yTie69v6vqxnkhKEClmZkcDHwFfAwPdPRc4HFgVPtZn/+vurYDdgNHAU2bWJ70hgZl1BToDB4TxZQNPAA+aWUtgBNAMeC/O2/cHJpe/cPdZwBzgnGTHXY/Uye89yRrHZ3Z3TSmcgG+AR2sosxvwOLAEWAr8E2gbs3w+8BvgHWADMB04OFz2N2AbUBIu+zqcfw0wC1gPLARuBTJi1tkKuB2YF5aZARxS1frixDwOuDHmdTbgwBkx86qMoZq4W4ZxfQusBsYCfWq5z38IFFf6vMeE22sF/AO4P853sCn8DBvC6dJw2e+Al9N9LNWFKeL3Xu13GB7PNxIk6A3ANGBfgiQ8BygCHgYya/p9AFcBX1SKsVf4XedH/C1s95nqy7EOLAOerTTvmHB+6535fjNpYOyWW+4CBqZwk1P85puvjVLQzPoCfYCf1lD0SWAr0D98/S+C/3hPiClzMXAKwUF4O8GPZA93v8rM9gHedvc/xJQvBI4j+DEOJDj45gMPhMsfAboAR4bz+wBezfqq+5zNYj7j7CgxVLOdh4FcYCiwBvgt8KqZDXD3rVHiAQYDM9x9WxhfPkGCvdfdN5jZIIJ9XsHdV5nZ6cB97t6r0vqmAZdG3HZC3HLLLRcCF+7CKkbffPPNo+Ot7+abbx6+C+utUM33HuU7vAA4mSAhPAa8QPAP0H4ECWEi8C7wFNX/Pp4Ebjezge4+JVx+ITDO3eeHr2v6LSTiM6fjWB8GzDazDu6+3MyaAncD17t7UW0/G9DwEgTBl1FXL9W0Dx+/q6qAmXUhyPp93X1NOO8XwCwz6+zuS8KiD7j7jHD5w8C1Zta6qgPB3Z+LefmFmT1BkAweMLMOwJnAPu7+bVjmm534fL81s18BOQQ/4EvdfWqUGOKtzMzaEfwX2dPdl4XzbgGuBYYQXKqLYjAwwMzWEhzz2cBLwK/C5W2AdXHeVwBMijN/HdA24rYTJZ9dO67HJXh9sar83mvxHT7o7jPD5U8B5wJD3X0jsNHMxgGDw8dqfx9m9hJwEXBNWL90AcEfW6D2x2FtP/PObCMRx7q7zzGz74C9gOXhe4sILoHtlIaYIKbUXCRt21sRPnYFZlZRpnv4+G3MvLkxy8oTxJKY5RvDxxyCA2IHZnYO8Atgd4LvvRnwabg4P3ycveM7a+WP7v4HM2tDcEZyRPgYJYZ4yv9znxr8zis05fv9FEUBcJW7PxTG0Rv4APglcBvBf2u5VbwvXny5BJcAUmk+8P4uvj+R64tV3fce9TuMPZ43AdvcfUWleTlE+308BvzLzP4HOBTIA54vL7wTx2E8dfVY/xrY08y+JkiKIzy81rQzGlyCiHq5Jx3cfbaZlVdwvl1FsUXhYz7B6TYEB1nsspqUxb4ws+4Ep+GnAa+7+xYzu53gDyB8/8djD+CrmtZXE3dfY2aXAnPN7BR3fylCDPG2s6A8rkp/LCILKw7bAONj4ptrZt/w/X79gu8vV8Q6APh7nPn7hO9JmfDy0Oi6uj6I/72TgO+wkii/jzeBzcCJwKnAv929GCL9FmqlLh3roVkEZxDDgOfcfcIurEutmNLgCuBcM/tTeDkJM+tgZjeY2VnuvpjgAL/DzPLC/1DuIDjQllSz3lhLCeoQyrUi+K5XAFvNbCjw4/KF7r4ceBa4z8zyLdAnplVG5fXVyN1XA3cSNAltUlMM8bYTxvVUGFdXgHCfnGpmrcLXvzOz+dWEMpjgrOqrsHxzMzsfODj8zAAvEpz+V9YesDjzjw7fI5VU/t6jfIe1XH+Nvw93LyOoxL6a4I/0ozGriHIc1jamunKsQ3AGcSpwEnDDrnwuUIJIOXd/i6B1UH9gmpmtBz4GOvD9Kf95BK0fZoXTWuD8Wmzmr0CBma01sxnhtd2bCa67rwWuB8ZUes/FBJfL3g+3/RLQKd76ahHH3QTNS8+PGEO87VxGcNCPC/fVNGAkQasRgB7seH091mCCSxPrzGwdwdnSucDx4XcB8AZQambDK733VmCMBe3djwAws34EZ1pPRdkBjVTF9x6+ruk7rK0ov4/HCOpYvo39Lzricbgz6sKxDsH+6AHcEiadXWK7cHlKJO3MbDZwpLtHvfxW1XqOBX7j7ofVUG4M8I67P7wr2xOprSjHupkNAd4iaPZbusvbVIIQEWkYzGwU8GN3PzQR69MlJhGRhmN/EtiSU2cQIiISl84gREQkLiUIERGJSwlCRETiUoIQEZG4lCBERCQuJQgREYlLCUJEROJSghARkbjqfXff7dq18/z8/HSHISJSr0yaNGmlu7evrky9TxD5+flMnDgx3WGIiNQrZragpjK6xCQiInEpQYiISFxKECIiEpcShIiIxKUEISIicaUsQZhZdzN7z8xmmtkMM7smThkzs3vMbI6ZTTWzQamKT0REtpfKM4hS4JfuvhcwFLjSzPpXKnMcwYDwewCjgH8kM6A5q1ezatOmZG5CRKTeSlmCcPcl7j45fL4emAl0rVTsFOBxD3wK5JlZ50THsnbzZo745z/Z4957eWjy5ESvXkSkQUhLHYSZ5ROMnfpZpUVdgUUxrwvZMYlgZqPMbKKZTVyxYkWtt9+6eXOWb9wIwMOTJ1OmYVdFRHaQ8gRhZq2A54Br3X1d5cVx3rLDX293f9DdC9y9oH37au8UryoGLhsUVG/MXbOGcfPn13odIiINXUoThJk1JUgOT7r783GKFALdY153AxYnI5bz9t2XZhkZQHAWISIi20tlKyYDHgFmuvudVRR7GTg/bM00FChy9yXJiGe3li05fa+9AHhu5kxVVouIVJLKM4hhwI+BI8xsSjgdb2aXm9nlYZn/AvOAOcBDwBXJDKj8MtOWbdt4YurUZG5KRKTeSVlvru7+EfHrGGLLOHBlaiKC4fn59GnbljmrV/Pw5MlcM2QIwYmOiIg06jupzYxL9t8fgBkrVvBpYWGaIxIRqTsadYIAuHDgQDKbBLtB90SIiHyv0SeITq1acVLfvgA8PWMG60pK0hyRiEjd0OgTBHxfWb1p61bGTJuW5mhEROoGJQhgRO/edM/NBXSZSUSknBIEkNGkCReHldWTlizhiyVJufVCRKReUYIIXbz//hVtcHVntYiIEkSFHq1bc2yfPgA8OW0am7ZuTXNEIiLppQQR49KwsrqopIRnZsxIczQiIumlBBHjpL596ZidDcCDuswkIo2cEkSMphkZXDRwIACfLFqkymoRadSUICq5vKCAJmF/TPdOmJDmaERE0kcJopKeeXmc0q8fAE9Nm8ZKdQMuIo2UEkQcVw8ZAkDJtm1q8ioijZYSRByH9+zJPh06AHDf559TWlaW5ohERFJPCSIOM+NnBx4IwKJ163hp1qw0RyQiknpKEFU4d8AA2rRoAaiyWkQaJyWIKmQ3a1YxmND7CxYwddmyNEckIpJaShDVuGLw4Ir+me797LO0xiIikmpKENXo1aYNJ4dNXp+cNo1VavIqIo2IEkQNyiuri0tLeeSLL9IcjYhI6ihB1OCIXr3o3749AH9Xk1cRaURSliDM7FEzW25m06tY3trMXjGzL81shpldlKrYqhPb5HVhURGvfP11miMSEUmNVJ5BjAaOrWb5lcBX7r4fMBy4w8yapSCuGp237760bt4cUJNXEWk8UpYg3P0DYHV1RYAcMzOgVVi2NBWx1aRVTJPX9+bPZ8rSpWmOSEQk+epSHcTfgL2AxcA04Bp3j3vB38xGmdlEM5u4YsWKlAT3syFDyAh7ef3zxx+nZJsiIulUlxLEMcAUoAswEPibmeXGK+juD7p7gbsXtA8rkJMtPy+PM/feG4D/zJjBt2vWpGS7IiLpUpcSxEXA8x6YA3wL7JnmmLbz62HDANjmzp3jx6c5GhGR5KpLCWIhcCSAmXUE+gHz0hpRJQM7dWJE794APPLFF6zYuDHNEYmIJE8qm7mOAcYD/cys0MwuMbPLzezysMj/Ageb2TTgHeA6d1+Zqviiui48iyguLeXvn3+e5mhERJInM1Ubcvdzali+GBiRonB22g/y8yno0oWJixdz74QJ/M/BB5PdrE60xhURSai6dImpXjAzfn3wwQCsLi5W9xsi0mBFShBm1s/M9oh5fYSZjTaz/zGzRpdkTttrL3q3aQPAHePHs3XbtjRHJCKSeFH/uI8GBgOYWVfgFYLmqD8Hfp+UyOqwjCZN+FV4FrGwqIj/zJiR5ohERBIvaoLYE5gYPh8JTHT3EcCPgXOTEVhdd8F++9EhOxuAP3/yCe6e5ohERBIraoJoCmwJnx8JvBY+nw10SnRQ9UFW06ZcHXbiN3XZMt6YOzfNEYmIJFbUBDEdGGVmBxEkiLHh/K7AqmQEVh9cMXgw2U2bAnCbut8QkQYmaoK4HrgC+Ah41t2nhvNPBhpt96ZtsrIYdcABAIybP59PCwvTHJGISOJEShDuPg5oB3Ry9/NjFj1CkDgarZ8PHUrTJsFuvOX999McjYhI4kRuourupcBWMzugfJwGd5/r7o267+vurVtXdAU+ds4cnUWISIMR9T6IVmFXGasJLil1C+f/w8xuSmJ89cJvDj2UZhkZAPxu3Lj0BiMikiBRzyD+D8gHDgSKY+b/Fzg9wTHVO91bt+bS8CzijblzGb9oUZojEhHZdVETxCkEA/hMJBj5rdxXwO4Jj6oeuiH2LEJ1ESLSAERNELsB8XpWbcX2CaPR6paby2WDBgHw5ty5fKKzCBGp56ImiInACTGvy5PCZQRdeAtwwyGH0Fx1ESLSQERNEL8BbjWzfxB0EX6Nmb1FMArcjckKrr7pmptbcV/EW/Pm8fHChWmOSERk50W9D+Ij4BAgF1gAHE9wB/XBYb2EhK6POYu4WWcRIlKP1eY+iCnufq677+nufd39bHf/MpnB1UddcnL4SXgW8c633/LhggVpjkhEZOdEvQ8it7op2UHWN9cfcggtMoPB+tSiSUTqq6hnEGuBNdVMEqNzTg6Xh2cR7377Le/Pn5/egEREdkLUBHE0wXjR5dMJBJXTC4ELkhNa/fbrYcMqziJueOcdjRchIvVO1ErqdypNY939VuDXwDnJDbF+6pyTw7VDhgAwvrCQF2bNSnNEIiK1s6vjSU8GhicgjgbpukMOoW1WFgDXv/22xq4WkXplpxOEmWUBVwORui81s0fNbLmZTa+mzHAzm2JmM8ys3tfu5rVowU2HHQbAN6tX8/DkyWmOSEQkuqitmNaY2eqYaQ2wHrgYuC7itkYDx1azjTzgPuBkd9+bYOzreu+nBQX0yssDghZN60tK0hyRiEg0mRHL/arS6zJgBfCpu8fro2kH7v6BmeVXU+RHwPPuvjAsvzxibHVa88xM/njEEfzo+edZvnEjd4wfz++GD093WCIiNYpaSf1Ipekxd381anKIqC/QxszGmdkkMzu/qoJmNsrMJprZxBUrViQwhOQ4a599OKBzZwBu/+QTlqxfn+aIRERqVmWCqOnmuCTcKJcJHEDQhPYY4CYz6xuvoLs/6O4F7l7Qvn37BG0+eZqY8ZejjwZg49atGppUROqF6s4garo5LtE3yhUCY919Y3hm8gGwX4LWnXY/6NWL4/r0AeDhyZOZtTKRJ18iIolXXR3E0SmLIvAS8DczywSaAUOAv6Y4hqS67aijGDtnDtvcueGdd3jhrLPSHZKISJWqTBDu/k4iNxSOaT0caGdmhcDNQNNwW/e7+0wzGwtMJagEf9jdq2wSWx8N6NiRCwcO5LEpU3hx1iw+XriQYT16pDssEZG4rLZdQJhZO4L/8Cu4++JEBlUbBQUFPnFi/elxvHDdOva49142l5YyuEsXPr30UpqYpTssEWlkzGySuxdUV6Y2vbk+YmYbgGXAokqTRNQtN5dfHnQQAJ8vXszoKVPSHJGISHxR76T+MzAYOAvYDPwYuAH4juD+BamFGw45hG65QeOv699+mzXFxWmOSERkR1ETxAnAVe7+GrANmODufyZIEhcnK7iGKrtZM+4YMQKAFZs2aeQ5EamToiaINgRDjQKsA9qGzz8mGIpUamlk//78ID8fgL9//jlTly1Lb0AiIpVETRDzgJ7h81nAmeHzU4DViQ6qMTAz7j3uODLMKHPnZ6+/rjEjRKROiZogHgcGhc//D7jSzEqAO4HbkxFYY7B3hw787MADAfhgwQKenjEjzRGJiHwval9Mt7v7XeHzt4H+BBXVBe5+dxLja/B+N3w4HbKzAfjlm2+yYcuWNEckIhKoNkGY2SVmll15vrt/6+7/cfcvkhda49C6RQtuO+ooABavX88fP/ggzRGJiARqOoP4G7DEzB4ys6GpCKgxOn+//RjarRsAd4wfz+xVq9IckYhIzQmiM/Abgl5WPzGz6WZ2bXg3tSRIk7DC2oCtZWVcrQprEakDqk0Q7r7W3f/m7oOAAuB94Cag0Mz+Y2bHpCLIxqCgSxcuGxS0A3hj7lzGTG9Q3VCJSD0UeUxqd5/s7lcCXQhujmsLvGZm85MUW6Nz61FHVVRYXzN2LCs3bUpzRCLSmEVOEOXcvQT4AphCcNNcp0QH1Vi1zcri3uOOA2Dlpk38/I030hyRiDRmkROEmbUys0vN7BNgOnA8cCug/qoTaGT//pzcrx8A/5o6lbFz5qQ5IhFprGpMEGZ2qJk9BiwB7gJmA4e7e393/4u7L092kI2JmfH3448np1nQo/pPXn1V90aISFrUdB/EbGAcsDfwK6Czu1/o7h+lILZGq1tubsW9EQuLirjx3XfTHJGINEY1nUGMBQa6+4Hu/oC7r09FUAI/KSjgkHC0uXs++4xPCwvTHJGINDY1NXO92t2npSoY+V4TMx466SSaZWTgwKUvv8yWbdvSHZaINCK1bsUkqbNnu3bcdNhhAMxYsYLbPtKVPRFJHSWIOu7Xw4YxoEMHAP7w4YcaN0JEUkYJoo5rlpHBwyefTBMztmzbxnnPP09JaWm6wxKRRiBKM9emZva0mfVORUCyowO7duW3hx4KwLTly9WqSURSosYE4e5bgRHALvUeZ2aPmtlyM6u2kyEzG2xm28zsjF3ZXkNz02GHUdClCxD0+Dpu/vz0BiQiDV7US0zPA6ft4rZGA8dWV8DMMoDbAPUxUUnTjAz+deqpZGVm4sAFL75I0ebN6Q5LRBqwqAliIXCjmb1kZjeZ2S9ipygrcPcPqHn86p8BzwG6OzuOfu3a8ZejjwaCG+h+9vrraY5IRBqyzIjlLgTWAPuGUywnGJt6l5hZV+BU4AhgcA1lRwGjAHr0aFxdQV0xeDCvzJ7NG3Pn8sTUqZzUty8j99473WGJSAMUdUzqXtVMuycolruA69y9xrvB3P1Bdy9w94L27dsnaPP1g5nx6Cmn0DYrC4DLX3uNxet1g7uIJF6tm7maWUczS0bz2ALg3+H4EmcA95nZD5OwnXqvS04OD554IgCri4u56KWXNAKdiCRcpD/0YVPXP5vZeuA7ID+cf5uZXZGIQMKzkXx3zweeBa5w9xcTse6G6PT+/Tl/v/0AeHPuXO4cPz7NEYlIQxP1TOBm4CTgPKAkZv4EgvqJGpnZGGA80M/MCs3sEjO73Mwur0W8EuOeY48lPy8PgOvfeYfxixalOSIRaUgsyqUJM5sLXOzu74dnEfu5+zwz6wd85u55yQ60KgUFBT5x4sR0bT7tPv/uO4Y9+ihby8rolpvLFz/5Ce1atkx3WCJSx5nZJHcvqK5M1DOILsCCOPMzid4SSpJgcNeu3HnMMQAUrlvH+S+8QJnqI0QkAaImiBnAYXHmnwlMSlw4sjOuHDyYkf37A/D6nDn8+eOP0xyRiDQEUf/7vwX4l5l1BzKAkWa2J/Aj4IRkBSfRWDh2xBdLlzJn9Wp+++67HNStG4fn56c7NBGpx6LeB/EKwdnCCKCMoNJ6D+Akd387eeFJVK1btOCZkSNpnpFBmTvnPPccyzZsSHdYIlKPRb6fwd3fcPfD3b2Vu7d090Pc/c1kBie1M7BTJ+457jgAlmzYwHkvvMC2srI0RyUi9VXU+yDmmdlucebnmdm8xIclO+uyQYP40YABALw9bx43jxuX3oBEpN6KegaRT1D3UFlzoGvCopFdZmY8cOKJ7NmuHQB//PBDnp5ebQ/rIiJxVVtJbWaxXXyfYGZFMa8zgCOB+UmIS3ZBq2bNePGssxjy8MMUlZRw0UsvscduuzGoc+d0hyYi9Ui1N8qZWfkFbAes0uKtBMnhl+7+alKii6Cx3yhXnTfmzOH4p56izJ1uublMvOwyOrZqle6wRKQO2OUb5dy9ibs3IRgPokP563Bq7u790pkcpHrH9OlTMX5E4bp1nPaf/2g8axGJrDbdfa9MdjCSeD8fOpQLwk79Plm0iJ++9pp6fhWRSCI3czWzK8xshpltMrPdw3nXm9mZyQtPdpWZcf+JJzK0WzcAHpsyhXs++yzNUYlIfRC1meu1wI3Ag2xfF/EdcFUS4pIEapGZyfNnnknXnBwAfvHmm7w5d26aoxKRui7qGcTlwGXufjcQexF7MqDxLuuBzjk5vHj22bTIzKTMnZHPPMOXS5emOywRqcOiJoieQLzG9FuBrMSFI8lU0KULo085BYB1JSUc9+STLFi7Ns1RiUhdFTVBzAMGxZl/PPBV4sKRZDtrn324Y8QIIOiO47gnn2R1cXGaoxKRuihqgrgd+JuZnUtQB3GQmd0M/BH4S7KCk+T4xUEH8fOhQwGYuXIlJ48ZQ/HWrWmOSkTqmqjNXB8Dfgf8CWgJPAFcClzt7k8nLTpJmttHjODMvYPqo48XLVLHfiKyg9r05vqQu/cEOgCd3L27uz+SvNAkmZqY8fgPf8jwcMyI52fO5JqxY3WPhIhUiJwgyrn7SndfnoxgJLWaZ2bywllnsU+HDgD8/fPPuU2j0YlIKOp9EG3M7G4zm2pmS81seeyU7CAlefJatOD1c8+lW24uADe88w5/nzAhzVGJSF0QdcjRxwnud/gnsIyg8z5pILrl5jL23HM5bPRoVhcXc9Xrr9MiM5NLBsVruCYijUXUBDEcONzdJ+/shszsUeBEYLm77xNn+bnAdeHLDcBP3f3Lnd2e1M7eHTrw5nnnceTjj1NUUsJlr7xC88xMztt333SHJiJpErUOYm4tylZlNHBsNcu/JUhC+wL/S9Cth6TQAV26MPa882jVrBkOXPDiizz7lW5zEWmsov7Rvwa41cz2M7N4I8vVyN0/AFZXs/wTd18TvvwU6LYz25FdM7RbN1770Y/ICrvkOOe553jl66/THZaIpEHUBDGHoEuNycAWM9sWOyUhrkuA16taaGajzGyimU1csWJFEjbfuB3Wsycvn3MOzTMyKC0r44xnnlHnfiKNULUjylUUMvsAaAPcT5xKand/LtLGzPKBV+PVQcSU+QFwH3CIu6+qaZ0aUS55Xps9m1OffpqtZWW0CJvEHtunT7rDEpEEiDKiXNRK6gLgQHeP12FfwpjZvsDDwHFRkoMk1wl9+/LvM87gzGeeYXNpKSePGcPTZ5zBqXvtle7QRCQFol5i+grITWYgZtYDeB74sbvPTua2JLrT9tqLZ888k2YZGWwtK2PkM8/w5NSp6Q5LRFIgaoK4EbjTzI4ys45m1jZ2irICMxsDjAf6mVmhmV1iZpeb2eVhkf8H7AbcZ2ZTzEzXjeqIH+65Jy+ffTZZmZlsc+fHL7zAg5MmpTssEUmyqHUQsb24xb7BAHf3nWrZlAiqg0idDxYs4ISnnmLDli0A3DliBD8/6KA0RyUiOyORdRA/SEA8Us8d1rMn75x/Psf+61+s2byZX7z5Jhu3buW3hx6KmdW8AhGpVyIlCHd/P9mBSP1wYNeujLvwQo5+4gmWb9zITe+9x/KNG/nrMceQ0WRX76UUkbqkVr9oM+tiZkPN7LDYKVnBSd20b8eOvH/hhXTNyQHg3gkTGPnMMxp0SKSBidqbaxczGwcUAh8D44D3YiZpZPZs147xl1xS0VX4C7NmccTjj7Ny06Y0RyYiiRL1DOIuYBvQH9gEHAqMBGZSff9K0oB1b92aDy+6iB+Egw59WljIwY88wtzVVfaoIiL1SNQEcThwnbvPImjFtMLdnyfoffV/kxWc1H15LVow9rzzOHfAAAC+Wb2agx55hAnffZfmyERkV0VNEFnAyvD5aoJhRyG4gU79QTdyzTIyeOLUU7nhkEMAWLFpE8NHj+b5mTPTHJmI7IqoCWIWsGf4fApwuZn1BK4E9K+iYGb86cgjuf+EE2hiRnFpKaf/5z/c/N57lGmca5F6KWqCuBvoFD7/PTACmAdcAfwmCXFJPfWTggJeOecccps3B+D3H3zAaU8/zfqSkjRHJiK1FelO6h3eZNaS4IxiobuvrKl8MulO6rpp1sqVnPLvfzN7VdDnYv/27Xnp7LPp0zZSzywikmRR7qSu8QzCzJqa2VIz27t8nrtvcvfJ6U4OUnft2a4dn116KcfvsQcAX61YweCHHtK4EiL1SI0Jwt23AlupNAaESE3yWrTg5bPP5vphwwBYu3kzxz35JLd99JHqJUTqgah1EPcCN5hZ1L6bRADIaNKEWz7r9UsAABYoSURBVI86ijGnn14xjOn177zDSWPG6KY6kTouam+urxDcC1EMTAc2xi5395OTEl0EqoOoP75YsoSRzzzD3DXB0ONdc3L49xlncEiPHmmOTKTxSUgdRGgl8BzwX2AhsKrSJFKj/Tt3ZtKoUYzs3x+A79avZ/jo0dz64Ye65CRSB+1UK6a6RGcQ9Y+7c//Eifz8jTco2bYNgGN69+aJU0+lfXZ2mqMTaRwSeQYhkjBmxk8HD+bTSy9lj7DZ6xtz57Lv/ffz2myNNitSV0ROEGZ2kZm9aWazzGxe7JTMAKXhGtipE5NGjeKcffYBYOmGDZw4Zgw/eeWVilHrRCR9onb3/T/AHcAkIB94kaCyui3waLKCk4Yvp3lznjztNEafcgo5zZoB8ODkyex3//18tHBhmqMTadyinkFcBoxy9xsI7on4W9hy6Q6gZ7KCk8bBzLhg4ECm/fSnFV2Hz1uzhsMee4zr3nqLktLS9AYo0khFTRDdgAnh82IgN3w+Bjg90UFJ49QzL4+3zz+fvx5zDC0yM3Hgz598QsFDD/FZYWG6wxNpdKImiKVAu/D5AuCg8HkfdIe1JFATM64dOpTJo0ZxQOfOAExfvpyDHnmEq19/XZ3+iaRQ1ATxLlB+M9wjwJ1m9h7wNPB8lBWY2aNmttzMplex3MzsHjObY2ZTzWxQxNikAdqrfXvGX3IJvx8+nGYZGTjB2Nd733cfr6qlk0hKRL2TugnQxN1Lw9dnAcOA2cADYX9NNa3jMGAD8Li77xNn+fHAz4DjgSHA3e4+pKb16j6Ihm/WypWMeuUVPoyptB7Zvz/3HHccnVq1SmNkIvVXlPsgIvWt5O5lQFnM66cJzh4ic/cPzCy/miKnECQPBz41szwz6+zuS6pb76ZNmxg9enRtQtnOwIEDGThwYMXrKVOmMGXKFAAuvPDC7cruynaqW191MeyMhviZ7h44kIn77sv/vPUWRSUlfPPVV9w4cya98vK47qqryGzy/clwfflMDfF70mcavdPbqW59yfhMUVSbIMJxH/4C/BBoCrwNXJ2kbr67AotiXheG83ZIEGY2ChgF0K1bNxYsWLDTG80PW82UW7t2bZXr25XtVLe+2sQQRUP9TJcNH86Jffty9dixrPzqK7q7U7pmDfs/8AB3H3ssR/TqtV0MiaLvKTp9pvr5mapS0xnELcCFwJPAZuAc4B/AyJ2OrGoWZ17c61/u/iDwIED//v29Z8+db2mbl5e3w+uq1rcr26lufbWJIYqG/Jk65+TwzMiRjB47ls8mTWJzaSnTly/nyMcf5/S99uL2ESPq3WeKsj59pprpM0VXeX1VqbYOwszmAr9193+Hrw8EPgZauPu22gYVXmJ6tYo6iAeAce4+Jnz9NTC8pktMqoNovDaXlvLX8eP5w4cfsmlrUA3WIjOT64YN49fDhtGyadM0RyhSdyWiL6buwIflL9x9AlAKdNn18HbwMnB+2JppKFBUU3KQxq1FZiY3HHooX191FT8aMAAIksYt779Pn3vu4cFJkygtK6thLSJSlZoSRAZQuVOcUiJWbscyszHAeKCfmRWa2SVmdrmZXR4W+S8wD5gDPARcUdttSOPULTeXJ087jY8uuohB4b0TSzZs4Cevvsre993H8zNnUt97LRZJh5ouMZUBbwGxdycdB7wPVAwHpgGDpK7YVlbGU9OmceN777GwqKhi/pCuXbntqKM4PGLlnEhDF+USU00J4rEoG3L3i2oZW8IoQUg8m0tL+cfnn/OHDz9kdXFxxfwRvXtz8+GHc3D37mmMTiT9djlB1AdKEFKdos2b+fPHH/PXTz+lOKbTv6N3352bDz+cYRruVBopJQiR0OL16/nThx/y0OTJbNn2fQO8I3v14ubDD+fQBDdPFKnrlCBEKilct47bPvqIhyZPrhjuFODwnj359bBhHNenD2bxbskRaViUIESqsHj9em776CMenDyZzTGXnvbp0IFfHXQQ5wwYQLOMjDRGKJJcShAiNViyfj23f/IJD06evN0wp11zcrhmyBBGHXAArVu0SGOEIsmhBCES0drNm3lg4kTu+uwzlm7YUDE/p1kzLthvP6488ED2bNeumjWI1C9KECK1VFJaypPTpnH7J58wc+X2fVIetfvuXDV4MCf27UtGk6hDqYjUTUoQIjupzJ3Xv/mGeydM4I25c7db1rN1a35aUMCFAwfSUeNRSD2lBCGSALNXreK+zz/nsSlTWBcz5Glmkyac1Lcvlw4axDG9e+usQuoVJQiRBNqwZQv/mjqVv02YwIwVK7Zb1jUnh4sGDuTi/fenV5s2aYpQJDolCJEkcHfGFxby8OTJPD1jRkVX4+UO79mTcwcM4Iz+/WmTlZWmKEWqpwQhkmTrSkp4evp0Hv7iCyZ89912y5plZHBi376cN2AAx++xB80za90JskjSKEGIpNC0ZcsYPWUKY6ZPZ0lMU1mAvBYtOG3PPTmjf3+O3H133YQnaacEIZIG28rKeG/+fJ6cNo3nvvqK9Vu2H1Ilr0ULTunXj5H9+3PU7rvrzELSQglCJM2Kt27lldmzeXLaNN6YM2e7/p8Acps356S+fTm5Xz+O7dOH3ObN0xSpNDZKECJ1yLqSEl6dPZtnv/qK1+fM2a4PKICmTZowPD+fk/v146S+fekZcWB5kZ2hBCFSR60vKeG/33zDszNnMnbOnO36gSq3b8eOHNenD8f07s2wHj1UbyEJpQQhUg+UlJYybv58Xv76a16ePZvCdet2KNOqWTOO6NWLY3v35pg+fdhd91rILlKCEKln3J0pS5fy0tdf8/qcOXz+3XfE+4X2ysvjiF69KqZO6vJDakkJQqSeW7lpE2/Pm8fYOXMYO2cOyzZujFuuf/v2HJGfz+H5+RzSo4cShtRICUKkASlzZ+qyZbw1dy7vzp/PBwsW7HAXd7k+bdtyaI8ewdSzJ73btNFIebKdOpcgzOxY4G4gA3jY3f+v0vIewD+BvLDM9e7+3+rWqQQhjdWWbduY8N13vPvtt7z77beMLyzcbrztWB2ysxnarRsHdevG0G7dKOjShVbNmqU4YqlL6lSCMLMMYDZwNFAIfA6c4+5fxZR5EPjC3f9hZv2B/7p7fnXrVYIQCRRv3crnixfz4YIFfLhwIZ8sWrTDTXrlmpgxoEMHhnTtSkGXLhR06cI+HTrQVC2lGo0oCSKVt3AeCMxx93kAZvZv4BTgq5gyDuSGz1sDi1MYn0i9ltW0KYf17MlhPXsCUFpWxtRly/ho4ULGFxbyaWEh89euBYLLVV8uW8aXy5bx4OTJADTPyGC/Tp04oHNnDujcmf07d6Z/+/a00J3ejVYqzyDOAI5190vD1z8Ghrj7VTFlOgNvAm2AbOAod58UZ12jgFEAPXr0OGDBggUp+AQi9d/SDRv4LEwW4wsLmbRkSdx7MMplmLFX+/bs17EjAzt1Yr+OHRnQsSMds7NVp1HP1bVLTCOBYyoliAPd/WcxZX4RxnSHmR0EPALs4+5lVa1Xl5hEdl6ZO7NXrWLi4sUV0xdLl1ZZ+V1ut6wsBnTsyD7t27NPhw7s06EDe7VvT1t1b15v1LVLTIVA95jX3djxEtIlwLEA7j7ezFoA7YDlKYlQpJFpYsae7dqxZ7t2nLfvvkBwaeqbVauYsnQpXy5bxpSlS5mydOl2TWxXFRczbv58xs2fv936OmRns2e7duwVrnOvdu3o164d3XNzNeJePZTKM4hMgkrqI4HvCCqpf+TuM2LKvA487e6jzWwv4B2gq1cTpM4gRFJj6YYNfLl0KTNWrGD68uVMX76cGStW1Hi2AcHYGL3btGGP3Xajb9u27LHbbvRp25bebdrQTckjLerUJSYAMzseuIugCeuj7v5HM/s9MNHdXw5bLj0EtCKosP61u79Z3TqVIETSp8yd+WvXMn35cmatXMnMlSuZuWIFM1eu3G787uo0bdKE/Lw8erdty+55eezepg35eXnk5+XRMy+P3bKyVN+RBHUuQSSDEoRI3ePuLN2wgZkrV/LNqlV8s3o1s8PHuatXs7WsymrFHWQ3bUrPMGH0yM2le+vWdM/NpUfr1nRv3ZquOTkaU2MnKEGISJ2zrayMBUVFzF29mnlr1jB3zRrmhdPcNWsin3nE6pCdTdecHLrm5gaPOTl0CafOOTl0btWK9tnZNNGZSIW6VkktIkJGkybs3qZN3B5p3Z3VxcUsKCpiwdq1zC+fioqYv3Yti4qKWLN58w7vW75xI8s3buSLpUur3q4ZHVu1onOrVnRq1YqO2dl0bNWKDtnZFc87ZmfTPjub3bKyVC+CEoSI1CFmxm4tW7Jby5YM6tw5bpkNW7ZQuG4dC4uKWFRUxKJ16/hu3Tq+W78+mNatY1Vx8Q7v2+bO4vXrWbx+fc1xAG2zsmifnU2H7Gzat2xJu5Yt2S0rK3iMeb1by5a0zcqidfPmDS6pKEGISL3Sqlmziqa5VdlcWsri9etZsn49SzZsqHhcumFDxeOyDRtYvnEj2+JcZneCpryriouZtXJlpLgMaN2iBW2zsmiblUWbFi1ok5VFXvPm5JU/b9GiYmrdvDmtY563bNq0zlXGK0GISIPTIjOzystYscrCS1rLNmxgWXiZakX546ZNwRQ+X7VpE6uLi+MmFAiSytrNm1m7eTPz1qypdcwZZuSGSSO3efOKKadZs2Aqfx7z2KpZM47p3TtplfRKECLSaDUxo114uWjvCOXL3CnavJlVxcWsjEkasdOacPma4mLWbt7Mms2bWVNcXGPLrW3uQdk4dSzVKbr+eiUIEZF0a2JGm6ws2mRl0adt28jvc3eKS0uDhFFcTFFJCUWbN1NUUsLazZsrnq8vKWHdli2sKympmIo2b2b9li2sLylhY5ybEpPZbbsShIhIkpkZLZs2pWXTpnTJydnp9ZS5s3HLloqEsWHLlqQ23VWCEBGpJ5qYBXUQzZvDLiSayNtL+hZERKReUoIQEZG4lCBERCQuJQgREYlLCUJEROJSghARkbiUIEREJK56Px6Ema0AFuzk29sB0XriSq26GhfU3dgUV+0ortppiHH1dPf21RWo9wliV5jZxJoGzEiHuhoX1N3YFFftKK7aaaxx6RKTiIjEpQQhIiJxNfYE8WC6A6hCXY0L6m5siqt2FFftNMq4GnUdhIiIVK2xn0GIiEgVGkWCMLNjzexrM5tjZtfHWd7czJ4Ol39mZvkpiKm7mb1nZjPNbIaZXROnzHAzKzKzKeH0/5IdV7jd+WY2LdzmxDjLzczuCffXVDMblIKY+sXshylmts7Mrq1UJmX7y8weNbPlZjY9Zl5bM3vLzL4JH+OOd2lmF4RlvjGzC1IQ11/MbFb4Xb1gZnlVvLfa7z0Jcf3OzL6L+b6Or+K91f5+kxDX0zExzTezKVW8Nyn7q6q/DWk5vty9QU9ABjAX2B1oBnwJ9K9U5grg/vD52cDTKYirMzAofJ4DzI4T13Dg1TTss/lAu2qWHw+8TjBO+1DgszR8p0sJ2nGnZX8BhwGDgOkx8/4MXB8+vx64Lc772gLzwsc24fM2SY5rBJAZPr8tXlxRvvckxPU74FcRvutqf7+JjqvS8juA/5fK/VXV34Z0HF+N4QziQGCOu89z9y3Av4FTKpU5Bfhn+PxZ4EizJA7TBLj7EnefHD5fD8wEuiZzmwl0CvC4Bz4F8syscwq3fyQw19139gbJXebuHwCrK82OPY7+CfwwzluPAd5y99XuvgZ4Czg2mXG5+5vuXhq+/BTolqjt7UpcEUX5/SYlrvBvwJnAmERtL2JMVf1tSPnx1RgSRFdgUczrQnb8Q1xRJvwhFQG7pSQ6ILyktT/wWZzFB5nZl2b2uplFGVc9ERx408wmmdmoOMuj7NNkOpuqf7Tp2F/lOrr7Egh+5ECHOGXSve8uJjj7i6em7z0ZrgovfT1axSWTdO6vQ4Fl7v5NFcuTvr8q/W1I+fHVGBJEvDOByk23opRJCjNrBTwHXOvu6yotnkxwGWU/4F7gxVTEBAxz90HAccCVZnZYpeXp3F/NgJOBZ+IsTtf+qo107rvfAqXAk1UUqel7T7R/AL2BgcASgss5laVtfwHnUP3ZQ1L3Vw1/G6p8W5x5O72/GkOCKAS6x7zuBiyuqoyZZQKt2bnT4Voxs6YEB8CT7v585eXuvs7dN4TP/ws0NbN2yY7L3ReHj8uBFwhO82NF2afJchww2d2XVV6Qrv0VY1n5pbbwcXmcMmnZd2Fl5YnAuR5erK4swveeUO6+zN23uXsZ8FAV20vX/soETgOerqpMMvdXFX8bUn58NYYE8Tmwh5n1Cv/7PBt4uVKZl4Hy2v4zgHer+hElSnh98xFgprvfWUWZTuV1IWZ2IMH3tSrJcWWbWU75c4IKzumVir0MnG+BoUBR+alvClT5X1069lclscfRBcBLccq8AYwwszbhJZUR4bykMbNjgeuAk919UxVlonzviY4rtt7q1Cq2F+X3mwxHAbPcvTDewmTur2r+NqT++Ep0DXxdnAha3cwmaA3x23De7wl+MAAtCC5ZzAEmALunIKZDCE79pgJTwul44HLg8rDMVcAMgpYbnwIHpyCu3cPtfRluu3x/xcZlwN/D/TkNKEjR99iS4A9+65h5adlfBElqCbCV4L+2Swjqrd4Bvgkf24ZlC4CHY957cXiszQEuSkFccwiuS5cfZ+Ut9roA/63ue09yXE+Ex89Ugj9+nSvHFb7e4febzLjC+aPLj6uYsinZX9X8bUj58aU7qUVEJK7GcIlJRER2ghKEiIjEpQQhIiJxKUGIiEhcShAiIhKXEoSIiMSlBCEiInEpQYhUw8xGm9mr6Y5DJB0y0x2ASLqYWU13if4TuIb4HaClhZndDuzj7gnrIlykKkoQ0pjF9gV0IkGHcbHzit29KLUh1Wgw8H66g5DGQZeYpNFy96XlE7C28jx3L6p8icnMxpnZP8zsDjNbbWYrzOwaC4at/buZrTWzhWb249hthR0b/trM5ppZcThU5XlRYzWzpma2hWAEtJvMzM1sRoJ2hUhcShAitXcusB4YAvwfcBfB2BOzCTpO+yfwsJl1iXnPHwg6qLuSYPjIW4EHzOyEiNvcBhwUPh9CcKZzyK59DJHqqbM+EcDMzgCecXerNH80wbjDJ4avxwHN3f2g8LUR9Ms/3t1PDuc1BTYCP3L3Z8PuoFcCI9z9w5h13wX0dffjI8Z4IvAUQW+2+uFK0qkOQqT2ppY/cXc3s+UE3VaXz9tqZmv4fkjI/gRdyo+tVDHelGDg+6j2B75UcpBUUYIQqb2tlV57FfPKL+GWP54ELKxhXdUZCHxRi/Iiu0QJQiT5vgJKCMbLfncX1rMf8HpiQhKpmRKESJK5+/rw/oXbwzqLD4BWwFCgzN0fjLiqTGDPsPJ7k7uvTU7EIgG1YhJJjZuA3wG/Ihii8i3gdOBbADO7MGy6ml/NOn5LMCZzIUErKJGkUismkTrAzG4BzgD2c/fSdMcjAjqDEKkrjgeuUnKQukRnECIiEpfOIEREJC4lCBERiUsJQkRE4lKCEBGRuJQgREQkLiUIERGJSwlCRETiUoIQEZG4/j9q333flbFMRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# How the contact rate and removal rate changes in time.\n",
    "\n",
    "A = 0.7\n",
    "B = 1.3\n",
    "alpha = 0.2\n",
    "mu = 1\n",
    "\n",
    "t_range = np.linspace(0, 20)\n",
    "contact_rate = rate(t_range, A, B, alpha)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(t_range, contact_rate, color = \"teal\", linewidth = 2.5, label = \"Contact Rate, $B(t)$\")\n",
    "ax.axhline(mu, linestyle = \"-.\", color = \"gray\", linewidth = 2.5, label = \"Removal Rate, $\\gamma$\")\n",
    "ax.legend(loc='lower left', bbox_to_anchor= (0.0, 1.01), ncol=2, \n",
    "            borderaxespad=0, frameon=False, fontsize = 13)\n",
    "ax.set_xlabel(\"Time, $t$\", fontsize = 14)\n",
    "ax.set_ylabel(\"Parameter Values\", fontsize = 14)\n",
    "#plt.savefig(\"imgs/exp_contact_rate.pdf\", bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:78: RuntimeWarning: invalid value encountered in double_scalars\n",
      "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed: 0.0\n",
      "Completed: 1.3333333333333333\n",
      "Completed: 2.6666666666666665\n",
      "Completed: 4.0\n",
      "Completed: 5.333333333333333\n",
      "Completed: 6.666666666666666\n",
      "Completed: 8.0\n",
      "Completed: 9.333333333333332\n",
      "Completed: 10.666666666666666\n",
      "Completed: 12.0\n",
      "Completed: 13.333333333333332\n",
      "Completed: 14.666666666666666\n",
      "Completed: 16.0\n",
      "Completed: 17.333333333333332\n",
      "Completed: 18.666666666666664\n",
      "Completed: 20.0\n",
      "130.63867902755737\n"
     ]
    }
   ],
   "source": [
    "# performing many simulations to obtain empirical probability of an epidemic.\n",
    "import time\n",
    "start = time.time()\n",
    "\n",
    "major_thresh = 30\n",
    "t_range = np.linspace(0, 20, 16)\n",
    "final_size = {}\n",
    "for t in t_range:\n",
    "    sims = []\n",
    "    for sim in range(5000):\n",
    "        t_hist, state_hist = time_dependent_gillespie_algo(0 + t, 50 + t, [4999, 1, 0], [A, B, alpha, mu])\n",
    "        if state_hist[-1][2] > major_thresh:\n",
    "            sims.append(1)\n",
    "        else:\n",
    "            sims.append(0)\n",
    "    final_size[t] = np.sum(sims) / len(sims)\n",
    "    print(\"Completed: \" + str(t))\n",
    "    \n",
    "end = time.time()\n",
    "print(end - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantenous rate\n",
    "\n",
    "t_range = np.linspace(0, 20)\n",
    "instant_rate = 1 - np.minimum(1, 1 / rate(t_range, A, B, alpha))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ODE approximation from MATLAB ode45\n",
    "\n",
    "ode_sol = pd.read_csv(\"data/exp_model_soln.csv\", header = None) # see MATLAB file labelled \n",
    "ode_time = pd.read_csv(\"data/time_exp_model_soln.csv\", header = None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAE5CAYAAABh4gz1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXhU5fXA8e+Z7CQhISyyQ9ghoqC4gFZUCoqKG1ZUXOtSq63WarWLVbG1rQrqz9ZWsbWuKK4tqBVR3I0LKiABFUhYAxhISEL2ZM7vj/cGhpCQZDKT9Xye5z7M3Llz3zcQ5sx73/ueI6qKMcYYcyC+lu6AMcaY1s+ChTHGmHpZsDDGGFMvCxbGGGPqZcHCGGNMvSxYGGOMqZcFC2OMMfWyYGGMMaZeFiyMMcbUy4KFMcaYelmwMMYYUy8LFsaYkBGRmSLyZpjO/biI/LEJ798tIoNC2aeOxIJFCxCR9SLywyae41IR+TCEfbpDRJ4O1fmag/fhUe59COSKyGIRGdGI9zf536GjEpFjReRjEcn3/u4/EpEjVPUZVZ3SCvr3rohcEbhPVRNUNbOl+tTWWbAwbd09qpoA9AG2AP9q4f60eyLSGXgV+CuQgvu7nwWUtWS/THhZsGhh1SMEEZktInkikiUiU2u8nikihd5rM0VkJPAwMN77Vr3LO/ZUEflKRApEZJOI3BFwnoEioiJyiYhsFJEdIvI777WTgd8CM7zzLff2J4nIv0Rkq4hsEZE/ikhEA/vdW0QWeN8614rIlQGv7XM5QUSOF5HNAc9v8dorFJFvRWRSfX+PqloCPA+MCTjPYBFZIiI7vZ/3GRFJ9l57CugPLPR+5pu9/Ud735h3ichyETm+If+OrY3MmtVbZs36gcya1TsMpx8GoKrPqmqVqpao6puquqLmiNf7nbtGRNZ4/55/8P5d0r3f0+dFJNo7dr/Rsvf+Ifv9fCJdRORVEcnxfv9eFZG+3mt3AT8A/ub92/6t5rm83+0nvfdvEJFbRcQX2I+6frc7KgsWrcNRwLdAN+Ae4F/ixAMPAlNVNRGYACxT1dXA1UC6N7RO9s5TBFwMJAOnAj8VkTNrtHUsMByYBNwmIiNV9Q3gT8B873yHesc+AVQCQ4CxwBQgcGhfa7+9154FNgO9gXOAPzXkQ19EhgM/A47wfuaTgPUNeF88cD6wNnA38GevDyOBfsAdAKp6EbARmOb9zPeISB/gNeCPuG/MNwEviUj3+tpvLWTWrATfrFnzgE3A+8Am36xZ82TWrIQQNvMdUCUiT4jIVBHpUs/xJwOHA0cDNwNzgZm4f4+Dcf9ujeUD/g0MwAX9EuBvAKr6O+AD4Gfev+3Pann/X4EkYBAwEff/5rKA1w/0u90hRbZ0B5rTgP/jnUa+Zc2G67mq5vs3XM8JAfvmAkNr7m+kDar6KICIPAH8HTgIKAT8wMEislFVtwJb6zqJqr4b8HSFiDyL+4/wn4D9s7xv4cu9EcShwOqa5xKRg4CpQLJ3fJGI3A9cBTxyoH6LSBQuKJ2mqqXAMhH5J3AR8HY9fxdVQAwwSkRyVHV9PcffJCI/AzoDG4AzAv4+1rI3eOSIyH3A7Qc414XA66r6uvd8sYgsBU7BBc5WT2CuiJz/q/HjOTE1lSVZWb7Z6enn44qcXRCKNlS1QESOBW4BHgV6isjrwJV1vOVuVS0AMkRkJfBm9dyBiPwP90WkUX+/qroTeKn6uTeaaND/b290PAMYq6qFQKGIzMH9flZfxqzr/+S2xvSzPbGRReuw5xdQVYu9hwmqWoT7pb4a2Coir8kBJnBF5CgReccbWud77+tWV1tAMVDXN84BQJTX7i5xl7oeAXrU12/cN/lc7z9itQ24a9sH5H3A/wI3AvheRJ4TkQNdSpntjawG4r5dDq9+QUR6eO/fIiIFwNPs//cRaADwo+qf1/uZjwV61dfv1kBmzeqtMOOm8eO5e/JkThoyhLsnT+bG8eNRmBHKS1KqulpVL1XVvrjRQW/ggToO3x7wuKSW540e9YhIJxF5xLuEVIAbRSV7gaA+3YBo3O9ktZq/n3X9bndYHWpk0YRv/nW+P3DkEQ6qughYJCJxuMsjj+Kux9ZWD3cebig+VVVLReQBDvzhuE9TNZ5vwk1YdlPVykZ2OxtIEZHEgIDRHzcBDe5yWaeA43vu0xHVecA8cROpjwB347711d151Y0icj3whIi86o2G/uz9XIeo6k7vktzfAt9W4zSbgKdUta5vyK3dYMB3YmrqPjsnpaZy78cf+3CXXLJD3aiqfiMijwM/ARY14VT7/F6ISM8DHHsj7ovBUaq6TUTGAF/hLj1C7f8/qu0AKnBfDlZ5+wJ/P00tbGTRionIQSJyunc9vgzYjbtMA+7bWd/qyUFPIu4bfamIHEnjLjtsBwZWT/J5l7zeBOaISGcR8XkTkxPrO5GqbgI+Bv4sIrEicghwOfCMd8gy4BQRSfE+EH4R8DMPF5ETRSQGKMV986yiAVR1Me7DsDqAJ+L+znZ58xG/quVnDrzv/mlgmoicJCIRXt+Pr544bQPWAf4lWVn77HzbPfcDIbltVERGiMiNARPK/XDzDp808dTLgTQRGSMisXjzS3VIxP1u7BKRFPa/vFjz33YPVa3C3Qxxl4gkisgA4Je4f39TBwsWrZsP9w0qG8jFzT9c4722BMgAtonIDm/fNcCdIlII3Ib7D9FQL3h/7hSRL73HF+OG66uAPOBFGn5J5nzcpaFs4BXgdu/DHOAp3AfDelxAmh/wvhjgL7hvf9twl71+24if417gZi/YzAIOA/JxE9cv1zj2z8Ct3iWnm7wgd4bXXg5upPEr2sj/E7399myB+bPT07l58WIWrV3LzYsXMyc9HYH5evvtoRpVFOImgD8VkSJckFiJ+10Nmqp+B9wJvAWsAQ60jugBIA73e/IJ8EaN1/8POMe7m+nBWt7/c9xIJtNrZx7wWFP6396J6oFGa8aYtkRmzUoQmKturssH+AXmK1ylt9++u6X7Z9ouCxbGtEPeZPYgIDOEIwrTgVmwMMYYU682cS3WGGNMy7JgYYwxpl4WLIwxxtTLgoUxxph6WbAwxhhTLwsWxhhj6mXBwrRK0owlMEXkYRH5fXO0ZUxbZcGihXgFVr4WkWIR2SYi/xCvMI/3+h0iUiGuYEyhiHwnIn8TkV4BxxwvIn7vgzVwG19P2+96aRBiwvkzNkW4SmBKLQV2VPVqVf1DqNsypj2xYNECRORGXCbVX+EKsByNy4C5uEZiwPleAaAU4CxcdtYvAgMGkO19sAZu6QdoeyB7s9aeHsIfq2Y7HSqjsTHtnQWLZual3Z4F/FxV31DVCq/Az7m4gHFhzfd4x2Tg8v3k0LSEbRfjEq89DlxSo2+Pe5dkFnujmfe8jJzVr6uIXCeuzOsOEblX9i1F+ZGI3C8iucAdXqbaW72aA9+LK2OZ5B0/wztPZ+/5VG+E1T2grSEB/fq7iPzPGzl9JCI9ReQBb4T0jYiMDejnr0VknfczrBKRs7z9dZWjrVnm9UpxpWBzxZWG7R3wmorI1eLKhOaJyEMiHbuCmukYLFg0vwlALDUyoKrqbuB/wOS63uilVv4vbmQQrItxqcKfAU4SVxEv0EzgD7g6GMvYm1a82lnAOFw21zOAHwe8dhQui2cP4C7gUm87AZenKIG9pS/nA+nAgyLSFVeh7ApVzamj3+cCt3r9KvPe+6X3/EXgvoBj1+H+jpJwgflpEel1gHK0e4jIibhstOfiMuxuAJ6rcdhpwBG4KoPn4kq/GtOudZhLBSdes7FZkmAt+Xv/+r5ldgN21FFQaCuuVvGBZOMuS1XrXf0NOUAfr8rePsSVwhwAPK+qO0RkHa7mxf0Bh72mqu97x/8OyBeRfl76bnAlMnOBXHHFlc4H/lndN1X9q/e4UkRmAvcFlND8DbBSRC7zfv5rgRXAu8BCVX31AD/3K6r6hXeeV4BrVPVJ7/l8XN1uAFT1hYD3zffaPRIXaOszE3hMVb8M6HOeiAwMKPH6F1Xdhaul8A4whv1TZBvTrtjIovntALrVcU2/l/f6gfTB1baolq2qyTW2/QKF5xJc/ePqNuZR41IUroYDsGe0k4srmbnf67hv3XW9hvdazdKVkbhaxngfuC/gynLOqaPP1RpcmlNELhaRZbK3NOrBNLxi4D599v4OdlJHyU0OXJrWmHajw4wsGvCNv7mk4y6jnE1AcSJx1fCmcoBCP978wDRccZhGEVeW9VwgQkSqP+xicHWLD1XV5d6+fgHvScCNYgJTXPfDFV0CV4oy8LWao7ds3EimWn+gEu+DXlwpzB8DzwIPAic39ueqyZtjeRSYhLvcVCUiy2hYuc39+uz9u3TFSm6aDs5GFs1MVfNx19H/KiIni0iUd4fSC8BmXBW5fXjHjMR9qPZk3+vzDXUmrjzpKNxlkzHASOAD3DxGtVNE5Fjvrqw/AJ8GXIIC+JWIdBFXSvN69q1yV9OzwA0ikuoFnj/h7vCqFFc282lccLwM6CMi1xzgXA0VjwsIOQAichluZFGttnK0geYBl4kr7Rnj9fnTgEtQxnRIFixagKreg/uQnA0UAJ/iLuFMUtWygENniMhuYBewAHc55HBVDfw231v2X2cxvZZmLwH+raobVXVb9YabcJ4ZcFlsHq6ecS5u/mRmjfP8F/gCN/n9Gm5iui6P4YLf+0AWrqb2z73X/gxsVtV/eD/zhcAfRWToAc5XL1VdhbuklY4LDKOBjwIOqa0cbeD73wZ+D7yEm0MaDJzXlD4Z0x4EVfzI+8bVG1cDN+cAd7CYNkREHsd9gN9ax+sKDFXVtc3aMWNMi2vwyEJEEkXkpyLyPpAPrMUVad8mIptE5FEROSJcHTXGGNNyGhQsROQGYD1uMnIx7v76McAwYDzuskUkbgXyG029lGCMMaZ1adBlKBF5AZilqivrOS4GuBwoV9V/HuhYY4wxbUdQcxbGGGM6lgNehvLy8XSrse/r6vw+xhhjOob65iwuBDrX2JeGW8y1DxFJEpGHQtUxY4wxrUd9wWLPqmcReV1E7sAteOpXy7GdgJ+ErmvGGGNai8ak+1gFHI8LIJ+JSCGwHPgKlwxuBG4RkzHGmHamvmDxJ7ykdap6E4CIlOFul+3N3rQRp3rnujlsPTXGGNNiGn03lIhE1pFeu1Xp1q2bDhw4sKW7YYwxbcoXX3yxQ1W719zfoMtQIpKqqlkA9QUKr2pY3xrJ55rdwIEDWbp0aUt2wRhj2hwR2VDb/oam+0gXkX+JyPgDNNBFRH6Km9s4I4g+GmOMaaUaOsE9Avgd8JqIVOGyjm7FZRHtgkt7PRL4DPiFqi4KQ1+NMca0kAaNLFR1l6r+Clct7KfAN0AykIorZvMEMFZVj7FAYYwx7U+jKuWpagnworcZY4zpIIIqfiQiMSJyqLfFNrUTXsW4b0VkrYj8upbXLxWRHK+u8jIRuaKpbRpjjGm4RgcLEZmMq+r2lbftFJF/ikiPYDogIhHAQ7j606OA80VkVC2HzlfVMd5mGW2NMaYZBTOy+BvwHjAE6A9cBAwFvhSRXkGc70hgrapmqmo58Bx2N5UxxrQqwQSL/sBvvA/3zar6sqpOBN7BBZLG6oMbqVTb7O2rabqIrBCRF0WkttxUIVdYBo99BVsKmqM1Y4xpvYIJFl8DfWvZ/0fgpCDOJ7Xsq7msfCEwUFUPAd7C3X21/4lErhKRpSKyNCen6WXBX1wNizLhqRVNPpUxxrRpwQSLecCjIjKkxv7uQF4Q59vMvlls+wLZgQeo6k5VLfOePgocXtuJVHWuqo5T1XHdu++3Wr3Rzh4JcZHw1Xb4aluTT2eMMW1WMMHiPmAwsEpEXhKRW0XkNtyH+I1BnO9zYKiIpIpINHAesCDwgBpzIacDq4Nop9GSYmD6SPf4qRVQ6W+OVo0xpvVp1DoLT3dgLHszzp4HDPdeu11EzsalLF+hqq/WdzJVrRSRnwGLgAjgMVXNEJE7gaWqugC4TkROxy0AzAUuDaLfQZk6BN7KhC2FsDjTPTfGmI4mJDW4RSQGGM3eADIGOFhVk5t88iCNGzdOQ5VIcGk23JsOCVHwwEmQuF+dQGOMaR9E5AtVHVdzfzAji/148wlLva3dObwXjO4BX38PL6yCH49t6R4ZY0zzCmoFd0cjAhcf4m7bWpwJm+xWWmNMB2PBooH6J8HkQeAHnlwOIbh6Z4wxbYYFi0b40SiIj4IV38MXVm3cGNOBWLBohM4xcI7dSmuM6YBCGixEpL+ItOsANGUw9EmErp1cOhBjjOkIQnI3VID1QIaIXKuq74f43K1CpA/umAiJ0W7i2xhjOoJQjwJ+DLwM3Bvi87YqnWMsUBhjOpaQBgtVfVxVb1fVo0J53tZqUwHc8xFsLWzpnhhjTHi16/mFcHt9DXyxDV5olkxVxhjTcho9ZyEidwGbVPXhGvuvBvqo6u9D1bnW7txREBMJZ41o6Z4YY0x4BTOyuAhXTrWmL4CLm9adtqVLHFx6qMtOa4wx7VkwwaIHUFtloZ3AQU3rTttV6bc0IMaY9iuYYLER+EEt+4/DFTLqcHYWw41vwl0fQGllS/fGGGNCL5hg8Qhwv4hcKSKDve0qYA4wN7TdaxtS4lwakLxSWPhdS/fGGGNCr9HBQlXn4ALGg8B3wBrg/3CV8tr1+oq6iMBFh7jHC7+D3JKW7Y8xxoRaULfOqupvgG7A0d7WXVV/raGopNRGjewOR/aGsip4PqOle2OMMaHVoFtnRWQBcKGqFniPazsGAFU9PXTda1suONhlo313A5wy1KU1N8aY9qChI4udgAY8PtDWYfVKdDUvFHjm65bujTHGhE6DRhaqelltj83+po+E9zbAsu2wYjsc0mFvJjbGtCeW7iPEOsfAmcPd46e/Bn+HncUxxrQnQaUoF5GewATcAr19Ao6q/j0E/WrTpg5xtbo35MMHG2HigJbukTHGNE0wuaEuBP4JCJDH3rkMvMcdPljERMKMNHhoKczPgKP7uH3GGNNWBXMZ6i7gHiBeVXuqaq+ArXeI+9dm/aA/pCbDzhJYsr6le2OMMU0TzPfdzsDjqmqJLQ5ABC4+BNbvcndIGWNMWxbMyOIZ4NRQd6Q9GtXdrbeItNsIjDFtXDAji18C/xGRScDXQEXgi6p6Zyg61t4UlEFFFXTt1NI9McaYxgsmWPwEOBnYAQxh/wluCxY1LNsGD3wKh/SAX45v6d4YY0zjBRMsfg/cqKr3h7oz7dWAJLfeoqzKjS6iIlq6R8YY0zjBBIsIoNb8UKZ2XeLg3slwUHxL98QYY4ITzNTrv4GZoeyEiJwsIt+KyFoR+fUBjjtHRFRExoWy/eZggcIY05YFM7LoBFwhIicBK9h/gvu6xpxMRCKAh4DJuEp7n4vIAlVdVeO4ROA64NMg+hyU7MJC1uXmMjglhd6JiSE5Z1YeLM6Cy8dAhN0lZYxpI4IJFiOBr7zHI2q8FkwmpCOBtaqaCSAizwFnAKtqHPcH3GLAm4Joo1F2l5dz1cKFzM/IQP0gPpiRlsbcadNIiI4O+rx+dRPd24pgcDJMsvUXxpg2otHBQlVPCHEf+gCbAp5vBo4KPEBExgL9VPVVEQl7sLhq4UJe+noNZ0WfjRQnknrcd8z5JB2AedOnB31en8C5afDgZ/DCajimP8RaGhBjTBvQGi6ESC379oxQRMQH3A/cWO+JRK4SkaUisjQnJyeozmQXFjI/I4PrJxxOVHEXduYK5/Y7lhvHj2d+RgbZhYVBnbfahL4wKNnV6359bZNOZYwxzSaoYCEiU0XkNRFZLSL9vH1XeAv1Gmsz0C/geV8gO+B5InAw8K6IrMeVcV1Q2yS3qs5V1XGqOq579+5BdAXW5ebiV2XSoFROPSYBgFc/3M2k1FT8qmTm5QV13moiMHO0e7zgW8gva9LpjDGmWTQ6WIjITOB54DtgIBDlvRQB3BxEHz4HhopIqohEA+cRcGuuquarajdVHaiqA4FPgNNVdWkQbdVrcEoKPhGWZGVx8vh4fD74aHkJb6zeiE+EQV26NLmNg3vAmIOgpBJeXh2CThtjTJgFM7K4GbhSVW8AApMJfgKMaezJvISEPwMWAauB51U1Q0TuFJFmr+fdOzGRGWlpzE5P557P32Fwqp8qP7zyQT4z0tJCdlfUzNHu+ttbmbC9KCSnNMaYsAlmenUokF7L/t24jLSNpqqvA6/X2HdbHcceH0wbjTF32jQA5qSnk1y6gUOZwtCqNB4+NXRVjPonwXH94b2NrubFdUeG7NTGGBNywYwssoFhtew/DljXtO60DgnR0cybPp1NN9zAK1dNoWuyj6qSWNasD22N1HPTIMoHH21y6y+MMaa1CiZYzAUeFJFjvOf9ROQS3BqIf4SsZ61A78REjhvYn2nHuktPr364O6Tn79YJpnhrLZ7NCOmpjTEmpBodLFT1HuBlYDEQD7wDPAw8rKoPhbZ7rcPU8fH4BD5cVkxuQVVIz33WCIiLhOXbYeX3IT21McaETFC3zqrq74BuuNXXRwPdVfX3oexYa9K9SyRHHRxHlR/e/CS0s9GJMXDGcPcPsSE/pKc2xpiQCXpRnqoWq+pSVf1MVUN7faYVOu1Yt+bitY924/eHdu5i6hC4bwqcOjSkpzXGmJBp0N1QIvJYQ0+oqj8Ovjut15GjYumeHMGWnEqWrSnjsOGxITt3bCT0Cs0ducYYExYNHVl0r7FNB87CVcobApwJnI27NNUuRUQIUye4POOvfhCegZQqfLYF0jeH5fTGGBO0Bo0sVHVa9WMR+Q1QAlymqkXevnjgX7ia3O3WKRMSePp/BXy4vJhdhVUkJ4a25N2K7TDnE0iOgbE9LcmgMab1CGbO4jrgjupAAeA9/gPw81B1rDXqkRLJkWmxVFbBGyGe6AY45CA4vCecPRIiW0OKR2OM8QTzkZQA9K5lfy9cYaR2bc9E94ehn+gWgZuPgZMGW7AwxrQuwXwkvQT8W0TOE5GB3nYe7jLUy6HtXutzVFoc3QImusOpPLRLOowxJmjBBIufAguBx3HpPTKBJ4DXgGtC1rNWKiJCOKUZJrpfWg3Xvg47isPShDHGNEowK7hLVPUaoCsw1ttSVPUaVe0QH22nTEhwK7qXF5NXGPqv/yKwuQAKyuHFmsVljTGmBTR0ncUC4EJVLfAe13YMAKra7GnFm1v1RPcnK0tZlF7EeVOCSrZ7QDPS4NMt8O4GOG0Y9A19E8YY02ANHVnsZG+p0531bB3CtB94yQXDsKIboGcCTEp1f+nzLcmgMaaFNXSdxWW1Pe7IjkxzK7qzcyr56rsyDh8RuhXd1c4eAe9tgM+y4budMKxryJswxpgGadINmiKSICIJoepMWxLhE04JqNEdDl3iXN4ogGdXuolvY4xpCUEFCxH5hYhsBPKBfBHZJCI3SPXERQdxyoTwpS6vdvowSIiCVTvcCm9jjGkJjQ4WInIPcAfwCDDZ2x4GbgPuDmXnWrvuXSI5erRLXf5GenhGF/HRLoU5wDwbXRhjWkgwI4srgCtU9S5VXeJtdwFXApeHtnutX+CK7qowTHSDW9GdEgvr8+FjSzJojGkBwc5ZrKhjX4dLUnHEqFh6do1g684qPl9VGpY2YiLhnFHu8fwMqPSHpRljjKlTMB/uTwLX1rL/p8BTTetO2xPhkz2jiwXvF4atneMHQO8E2F4ES7LC1owxxtQqmGARA1wqIt+IyOPethr4MRApIg9Wb6Htaus1dUICkRHwaUYp23ZWhqWNCB+cdzAMTILeVijJGNPMgqmYMAL40ns8wPtzm7eNDDiuw0zFdkmM4LixnViytJjXPtzN5Wckh6WdI3u7rWPdc2aMaQ0aHSxU9YRwdKStO/24BJYsLeb1j3dz8alJREWG/hO9ZpBQtcBhjGkeHW5COlxGD45hYK8o8gr9fLgsvPkUdxTDQ59bGhBjTPNpcLAQkY9FJDng+Z9FJCXgeTdvoV6HJCKcfpw30R2m1OXVdpXC+xth0TooDc8UiTHG7KMxI4ujgeiA59cCgRfnI4A+oehUWzX5yHhiY4Tla8pYv7UibO0MSYFLD4W//NDqdBtjmkdTLkPZ1fIa4uN8TD7CFUZa+EH4bqMFlzPqoPiwNmGMMXvYnEWIVV+KevOTIkrKwr96TtVlpDXGmHBqTLBQ9r8dtsPcHttQg/tGMyo1mqJS5e3PwzvR7Ve47V23bdgV1qaMMR1cY4KFAE+LyAKvWl4s8GjA8yeD7YSInCwi34rIWhH5dS2vXy0iX4vIMhH5UERGBdtWczhzols19593C9EwZv7ziZu/UOA5uzPKGBNGjQkWTwDZ7K2I9zSwKeB5NkEEDBGJAB4CpgKjgPNrCQbzVHW0qo4B7gHua2w7zWniYZ1I6ewjM7uC5WvKwtrWmSMgNgK+3Abf7AhrU8aYDqzB99KEsULekcBaVc0EEJHngDOAVQFtFwQcH08rv/wVFenyRT35egGvvFvImGGhr6JXLSnG1eh+cbUrkHTHRFuoZ4wJvdYwwd0HN0KptplabsEVkWtFZB1uZHFdM/UtaNN+kEiEDz5aXsL23PAuhjhtKCRGwzc7YZkVSDLGhEFrCBa1fQ/eb+Sgqg+p6mDgFuDWWk8kcpWILBWRpTk5OSHuZuN0TYpg4mGd8Cv89/3wLtKLi4KzRrjHVn7VGBMOrSFYbAb6BTzvi5v/qMtzwJm1vaCqc1V1nKqO6969ewi7GJyzT3AT3a9/tJuy8vDeRjt5EHSLgw358NGm+o83xpjGaA3B4nNgqIikikg0cB6wIPAAERka8PRUYE0z9i9oIwdGM7x/NAVFfpYsDe9ttNERMN3L+fv8KiuQZIwJrUYFCxGJEpFPRWR4qDqgqpXAz4BFwGrgeVXNEJE7ReR077CfiUiGiCwDfglcEqr2w0lEOOt4t0jv5TDfRgswMaBA0jvrw9qUMaaDaVSwUNUKIJUQ342kqq+r6jBVHezV80ZVb1PVBd7j61U1TVXHqOoJqtpmVhUcf3g8XRJ9rNtcwdfrwnsbbYQPZqS5xy+thjJLMmiMCZFgLkM9AVwZ6o60V9FRwqle2dV5i/P4YMMGsgvDlzfqqD4wKBnySuHNzLA1Y7l9ObQAACAASURBVIzpYILJWRoPzBSRycAXQFHgi6ra6m9rbW6Tjorh6Tf8fPp1GQ9sfZGKyGJmpKUxd9o0EqKj6z9BI4jABaNhdQ5MSg3pqY0xHVgwwWIke8uqDqrxmt20WYvffPA622NSOKh0MNf1n44MWMvs9HQA5k2fHvL2RvdwmzHGhIqVVQ2z7MJC5mdkcM3448l4B1ZmRDL/gkkoMCc9ndlTptA7MTFs7ZdVujuj4kM7gDHGdDCt4dbZdm1dbi5+VU4b25vDR8RSWqYsfL+QSamp+FXJzMsLW9vLt8H1iyzJoDGm6YKqsyYiB+Eq5Y3CXXpaBfxdVS3ZRA2DU1LwibAkK4sZk4/li29KeendQlJZj0+EQV26hK3tlE6uBGvWLqjyu7uljDEmGI0OFiJyDPAGsB1I93bPBG4QkZNUNb3ON3dAvRMTmZGWxuz0dPxHQ4/uQ/k+Bz55N4cZR6SF9RJUv87wxxNgcBdLLmiMaZpgRhazgWeBq1XVDyAiPuBhYA4wIXTdax/mTpsGwH2fpNOtbBtpTGRU5WE8fOrAsLc9JCXsTRhjOoBgLkyMAeZUBwoA7/F9wNhQdaw9SYiOZt706Wy64Qae/8lEunXxUVUcx4pvqpqtD9uL4LU2kSTFGNMaBRMs8nGruGtKBay45wH0TkxkYmp/ZkzqDMD8twrqeUdolFXCb9+GJ1fAt1YgyRgThGCCxXPAv0Rkppf8b6CIXAg8irs8ZepxyoQEEjv5yMgsZ2WYU4AAxETClMHu8bMZlsLcGNN4wQSLm4EXgceAtcA64J/AC8B+9bPN/uJifZwx0aUAmb+4eUYX04ZBQhSs3mEFkowxjdfoYKGq5ap6PdAFN38xFkhR1RtUtTzUHWyvzpyYSFQkfLSihI3bKsLeXqcoV68b4DkrkGSMaaSg77xX1WJV/VpVV6hqeIs1tEMpnSM4+Wg3unj+7eYZXUwZBCmxsD4f0jc3S5PGmHaiQcFCRBaISOeAxwfa5onIb0UkKbxdb/t+9MNERGDxp0XszA//nVExkXCOVyBpfoYVSDLGNFxDRxY72ZskcGc9G7gU5k+FrpvtU98eURx7aBwVlfDKu+FLWx7o+IHQKwG2FcF765ulSWNMO9CgRXmqelltj+siIqNw5VJNPWZM7swHy0r47/uFXHBSZzrFhjcnR3WBpAc+hRdWw7H93YjDGGMOJKhPJhGJFJEJInKeiFwcsF3kHfIttpK7QUalxjB6SAxFJcqrH+5uljaP7gOpXoGkReuapUljTBvX6GAhIiNwtbLfB57B3Tb7OG6dxUMAqlqlqstD18327fwp3iK9xQWUlIZ/IkEELjjYPf7Pt7Db7mEzxtQjmJHFA7gKeUlAMa4Y0jhgGRD6Sj4dwFFpsYxKjSav0M8r7zXP3MXoHpDWHYoqYNm2ZmnSGNOGBRMsjgD+qKpFgB+IVNUvcYv15oSycx2FiPDjackAzF9cyO6S5hldXHoo3HWCm7cwxpgDCSZYCG5EAZAD9PEebwaGhKJTHdHY4TGMGRpDYbGfF5tp3UX/JMtKa4xpmGCCxUrgUO/xZ8AtIjIRmIVL/2GCICJcNs0tTXlxSSH5u5svIy3A6hyXmdYYY2oTTLC4Cze6ALgV6Ae8A0wBrgtRvzqk0UNiOWJULMWlyvy3mmfuAlzq8jved2lAjDGmNsHkhlqkqi97jzNVdRTQDThIVd8Ncf86nMtOc6OLV94pJLcZVnUDHNUHkmKgb2fLGWWMqV1IVoCpaq6qfcyEwoiBMRxzSBxlFcq8N5tn7qJbJ/jbVJg+0sqvGmNqF97lwiYol3qji4UfFPJ9bmWztBkd0SzNGGPaqAYnehCRBQ05TlVPD747BmBw32hOOLwT73xRzDNvFHDDBc1zy1KlH5ZkwYrv4cajbZRhjNmrMSOL04DRNCyRoGmiS05Nwifw+se7yd7RPKOLSj+8uBo+z3abMcZUa0ywmA3EAMfhquP9XlUvq7mFpZcdUP+eUUw+Kp4qPzz1en6ztBkbCWdXF0jKgCpLYW6M8TQ4WKjqzbjbZG/ApfdYIyL/E5FzRCSqKZ0QkZNF5FsRWSsi+5VmFZFfisgqEVkhIm+LyICmtNdWXHRKEhE+V++iOarpAfxwEPToBFsK4b0NzdKkMaYNaNQEt5cgcIGqngmk4tZX/BHYIiIJwXRARCJwCQinAqOA870U54G+Asap6iG4+t/3BNNWW9O7WySnTEjAr/D4q80zuoj0UpgDvLAKyprnCpgxppVryt1Q8UAykADsZm9xpMY6EljrrdkoB54Dzgg8QFXfCSjd+gnQN8i22pwLp3YmKhLe/bKYVVllzdLmMf1gYBLkWgpzY4ynUcFCROJE5BIReR/4GhgAXKKqg7zEgsHoA2wKeL6ZvfmmanM58L8g22pzuneJ5NxJLoX5g/PzqPKHfzmLpTA3xtTU4GAhInOBbcDPgWeB3qo6U1XfbmIfartBs9ZPRBG5EDdfcm8dr18lIktFZGlOTk4Tu9V6XHByZ7onR/DdxnJe/2hvgaTswkI+2LCB7MLQpwY55KC9KcwXfBvy0xtj2pjGFNS8AtgIbMXNL0yVWm7ED2KdxWbcxHm1vsB+N26KyA+B3wETVbXW6zGqOheYCzBu3Lh2s6I8LsbHT8/pwp3/3MG/FuRz+MFR3LTkNeZnZOBXxSfCjLQ05k6bRkJ0dEjarB5d/O4d+N9aOHkIpMSF5NTGmDaoMcHiSYKflziQz4GhIpIKbAHOAy4IPEBExgKPACer6vdh6EOrN3FsHGOHx/DVt2Vc+feVvFGVwU3jx3NiaipLsrKYnZ4OwLzpoas/NSTFlWD9ZIub7P7J4SE7tTGmjWlwsFDVS8PRAVWtFJGfAYuACOAxVc0QkTuBpaq6AHfZKQF4wRvNbOxoK8VFhJ+fm8KVd22leEt3fnrC8dw9+TgAThoyBAXmpKcze8oUeicmhqzd89Lgsy3wzno4bSj06RyyUxtj2pBWkRtKVV9X1WGqOlhV7/L23eYFClT1h6p6kKqO8bYOFSiqDewVxZGHK4Kw65v++AMmuyelpuJXJTMvL6Rt9kqEE1OheyfIb56bsYwxrVCrCBam4S47rQvlvhKyt/pY/NneG9DezsrCJ8KgLl1C3ubM0XDfFBjVPeSnNsa0EY2ZszCtwJDunek+4kvyVw1izvzvKU3Yxkdbs5iTns6MtLSQXoKq1qlJ6/ONMe1Bg0YWInKciFhgaSWevHI8UUmFVJZFctPja/YEirnTpoW13cIyeHI5fLUtrM0YY1qhhl6GegdIARCRTBHpGr4umfokxsTw0DVDEYEBJWl8fN51zJs+PWS3zdblw03w2lp45murqGdMR9PQYJGHywUFMLAR7zNhMqRfNKcf5/JGPfdaGc1RqPCHqTCxP1w7zmpdGNPRNPTS0kvAeyKyFbfWYqmI1FogWlUHhapz5sAuOy2Jd78oZvmaMt76rJjJR8WHtb2oCLjmiLA2YYxppRoaLK4GFgBDgfuAfwOhzzFhGqVzfARXnZnMvU/n8uDzuRw6NIYeKc03tbSzGLp2arbmjDEtqEGfLOqucbwGICKHAnNU1YJFK3Dy+Hg+WFbMJytL+cuTO5l9XQ98vvBeIyqrhNnpsDYXHjwZEmPC2pwxphVo9NyDVw2vk4jcKSIvisgLInKHiBwUhv6ZeogIN13YleQEH8u+K+Old8Ifw2O8rxjFlfDKN2FvzhjTCjQ6WIjIMcBaXP6mEqAUuBBXOW98aLtnGiKlcwQ3XZgCwD//u4vMLeHPKV6dwnzROvg+2OT0xpg2I5i7mmbjUpQPU9WLVPUiYBiuaNGcUHbONNyEQzpx6jHxVFTCn/69k/KK8N4dldoFju0HlQrzM8LalDGmFQgmWIzBzVn4q3d4j+8DxoaqY6bxrpnehT7dI8nMruCxhbvC3t65aRApbv1FVmhTUhljWplggkU+e9dcBEoFwv8JZeoUF+vjN5d2xeeDF94uZNl3pWFt76B4OGmwe/y0LdQzpl0LJlg8B/xLRGaKSKqIDPQq2D2KuzxlWtCo1BguPLkzqvCXJ3ayu9hf/5ua4OyREB8FK3Ng+fawNmWMaUHBBIubgReBx3AT3euAfwIvAL8OXddMsC6cmsTwAdF8n1fFg/Nzw9pWQjScNcI9fvpraIYS4caYFhDMrbPlqno90AU3fzEWSFHVG1Q1/LfhmHpFRgi/vbQrsdHCW58X887S8N6udNJgV+9iUwG8tyGsTRljWkjQOZ5UtVhVv1bVFapaHMpOmabrd1AUV5+dDMB983LZsLUibG1FR7iKeuDujCqrDFtTxpgWYgkB27FpP0jguLFxFJUqv3s4h/zdtabzColj+kFqMuSVwutrw9aMMaaFWLBox0SEWy7uytB+UWTnVDLr0R1UVIZnUkEELjrEZaY9YWBYmjDGtCALFu1cXIyPP17dna5JESxbU8aD83PDls48rTtceRgkx4bl9MaYFmTBogPo3iWSO3/Sjego4bWPini5GfJHVfpdZT1jTPsQTG6o/4jIaSJigaYNGTkwhlsucvmj/vHSLj7NKAlbW1l58KvF8MiXYWvCGNPMgvnALwLmA5tF5E8iMjTEfTJhcsK4eC4+pTN+hT/8awdZ2fve6ZxdWMgHGzaQXdi0kUdyLOwsgU35UBy+m7CMMc1Igrl+LSKdgZnAZcA44EO8hXmqGr6vrI0wbtw4Xbp0aUt3o9Xx+5U/PLaT974splfXCP5+S08ioqu4auFC5mdk4FfFJ8KMtDTmTpsWdF3vNbnu7qhIG38a06aIyBeqOq7m/qDKqqlqAfAP4B8ikgZcATwC/FVEngMeUNXVTemwCQ+fT7jl4hS27qjku43l3DZ3B9/3+4jnV2dw0/jxnJiaypKsLGanpwMwb/r0oNoZmhLKXhtjWlqTvveJSG/gDOA0oBKXBqQfsEJEbmp690w4xEb7+MPV3eiaFMHXa8vI+LgHvzxqPHdPnsxJQ4Zw9+TJ3Dh+PPMzMpp8SSq/DJ5daQv1jGnrgpngjhKRc0TkdWADcCZwD9BLVS9X1VNwl6huDW1XTSh1T47kT9d0Jy5W6VbWn4KVQ/dZgzEpNRW/Kpl5Tcs9/sAn8J9vbaGeMW1dMCOLrbhLTuuAw1X1SFV9VFV3BxyzGLAKB63c0H7R/PaqzlRIGWvW+bjj0R17iia9nZWFT4RBXbo0qY2zR7o///MN5Ic3Y7oxJoyCCRY3AH1U9eequqK2A1Q1T1Vrq3lhWpljRnSh5xFrqZBS0r8u4dr/28BNb7zFnPR0ZqSl0TsxsUnnH90DDu8JpVUwf1WIOm2MaXbBBIv3gP2WW4nTv+ldMs3t3xdM2hMw1mX6WLwojnNHHszcadNCcv6Zo90v2pIs2JgfklMaY5pZMMEiC+hey/4U7zXTxiRER/PCpadwz3Xd6RQHXcv6MfD7Y4mWqJCcv09nmDwIFHjm65Cc0hjTzIIJFoL7f19TAhDUVWkROVlEvhWRtSKyXwElETlORL4UkUoROSeYNkz9xg/vwv/d0JOkBB+fZpRy29ycPXMYTTV9FMRFwrLtsHxbSE5pjGlGDQ4WIvKgiDyICxR/rn7ubQ/hbptd1tgOiEgE8BAwFRgFnC8io2octhG4FJjX2PObxhncN5o51/cgKcHHZxml/Pbv31MYgtKsSTFwtldR7ymrqGdMm9OYkcVobxNgZMDz0cAQ4EvcB3pjHQmsVdVMr9Lec7i1G3uo6npvMj28BaUNAIP6RHPfL3rQJdHHl9+Wce0929i4rfF5O2qmDzl5CPTwKuotsQuWxrQpDV7BraonAIjIv4HrvVXcodAH2BTwfDNwVIjObYKU2juah27uye8fyWHd5gquvWcbv/txN44+OK7e9+4uL68zfcgFo6N54FN4fhVM6AedQjMtYowJs2BqcF8WwkABbqSyXzNBnUjkKhFZKiJLc3Jymtgt07NrJA/eeNDeanv/yOG5xQX11sOoDhQ3jR/PGzNncpO3GvyqhQs5ug8M7+pWdr9sCWGMaTMaNLIQkQXAhapa4D2uk6qe3sg+bMalCKnWF8hu5Dmq254LzAWXSDCYc5h9xcX4uO3ybjz9RgGPv5rP3Fd2kbmlnBsvSCEmev/vGtmFhXsCxd2TJwNw0pAhKDAnPZ3ZU6ZwyaGJ/HYJbC4AVVdlzxjTujV0ZLGTvd/2c73ndW2N9TkwVERSRSQaOA84YEAyzcvnEy4+JYk7ruxGbIzw1mfF/OL+78nZtX/Cp3W5ufhVOTF13zWZgelDBneBv5wItxxjgcKYtqJBIwtVvSzg8aWh7ICqVorIz4BFQATwmKpmiMidwFJVXSAiRwCvAF2AaSIyS1XTQtkPU7/jxnaiT/dIbn04h283lHPN3du55eIUxo3cO48xOCUFnwhLsrI4aciQPftrpg9JbVoWEWNMM2tQPYv6Lj0FUFU9o/7Dws/qWYTPrsIq7nh0ByvWuoX8036QwNVnJRMX6waqF7z0EvMzMrhx/HgmpabydlbWnvQhNVOeb8yHxZlw2Rjw2SjDmBbX1HoWwVxeMu1UcmIEc67vwXOLC3jitXwWfrCbz1eVcPNFXRkzLHZPmpA56enc+/HH+9wNFcivMDsdthfBwCSYNKglfhpjTEMEVSmvLbCRRfNYt7mcu5/cydrNbh3GWccncMUZycTF+MguLCQzL49BXbrUmZDwk82wKgd+NAoSY5qz58aY2tQ1srBgYZqsskp5+n/5PPNGAVV+6NM9klsu7srBg+3T35i2pknBIsy3zoaFBYvm991GN8rIyq5ABM44LoFLTk0iKSGiQe+vqILiSpcaxBjTMuoKFsHcOnug22ZtbqMDG9Y/mn/c0pOZJ3VGgP+8t5sLb8vmuTcL6k1IuC4PblwMj3zRPH01xjSOXYYyYbFuczlz/7OLz1e5RMQ9UiK44vRkThzXCV8ttz3tKoVfLIKSSvj1MTC25/7nzC4sZF1uLoNTUppclMkYU7uwzFmISAJAjZKqrYIFi9bh81UlPPLKLjK3uAnwYf2jufrsZMYMi93v2Fe/cxlpeyXA7MkQ6Y17D5RrKiE6ujl/HGPavaZehqp5sl+IyEYgH8gXkU0icoOIrcc1+zpiVByP/KYnv7ooha5JEXy3sZxfPvA9v/3792Rk7ltw8eQh0DsBtu6G/63du/9AuaaMMc2jwVlnq4nIPcBVwL1Aurd7PHAb0Au4OWS9M+1ChE+YOj6B4w/rxItLCnnuzQI+WVnKJytLOXhwDOdOSmTCIXFE+oRLDoU/fwQvroJj+0FJZf25puySlDHh1+hgAVwBXKGqLwbsWyIi3wKPYMHC1CEuxsdFU5M49ZgEXn6nkAXvF7JyXRkr15XRt0ckP5qUyJSj4jm8l48vtsKzK2F0j7pzTd378cdk5uVZsDCmGQR1GQpYUce+YM9nOpCUzhFccUYy8+/qw7XnJHNQSgSbv6/k/mfzOP/WbBK25yPlVby3EXzSbU+uqUA1c00ZY8IrmJHFk8C1wPU19v8UeKrJPTIdRlysj+kndubMiYm8/1Ux898q5LuN5by4OB/x5ROVEsdTRfH8aMTBzE5PR2G/XFM2qjCmeTR0Ud6DAU8jgQtxNSc+8fYdBfQGnlHVa0LdyWDY3VBtj6qyfE0ZLy0pJH1lCX6viG5crCDdtvFByacUROzE57O7oYwJl6YmEhxd43n10qkB3p/bvG1EcN0zBkSEMcNiGTMsltyCKv75ZhGLPyuiZHcFbD6IwzmdHt3hxMPiOPXIZBKirSarMc3FFuWZVksVbntXydhQwYCy3WStLaagyL/n9YG9ojjmkDgmHBrH8P7RtS72M8Y0jiUSNG3S+l2QkQMnDQb1K5+sLOG9L4v5dGUJRaV7f3e7JkUw4ZA4jjkkjkOHxtRa8jWQrQY3pnYhDRYiEgkcCfQH9rlorKpPBtvJULJg0b5VVCrL15Ty0YoSPl5eQs6uqj2vRUVCWmoMY4fHMnZ4LMMHRBMV6UYd4VoNbsHHtBchCxYiMgJYCKQCAlTh5j4qgDJV7dz07jadBYv2Z1cpbC2Ekd333a+qrNlUwUfLi0n/umRPbY1qsTHC6MExjB0Wy4ubP+Y/m5Zy44SjODE1lSVZWcyuo4pfQ1gqEtPeNHWCO9ADuAnuMbhJ7TFAEvAP4NamdNKYumwpgN+9A9ERcP8UiA/4HBYRhvWPZlj/aC6blkz+7iqWrynjq+9KWfZtKRu2VfL5qlIvqeEoJvpGUPhVLN/lRTN5YC9KRkfx0Ir3g1oNHpiKJDD4AEEFH2Naq2CCxRHARFUtEhE/EKmqX4rIzcBfgUNC2kNjgN6J0D8JOkVCeRXEH+DYpIQIjhvbiePGdgJgZ34Vy74r5X9f7OCDjHziq5L5Zn0536wv5xV2A4OYIL2545E8jhhSxaA+UQzqE0W/g6KIjKh70jy7MHypSOyylmltggkWAhR7j3OAPsC3wGZgSIj6Zcw+ROA3x0BspHvcGF2TIph0RDwjR/i5//4n+MW4Y7hg0DGszirjmw3lfPZtIZTFsjYT1mYW7HlfVCT07xnFoN5RDOoTzcDeUfTtEUnPlEgiIoR1uaFPRWJzKqa1CiZYrAQOBTKBz4BbRKQKuBJYe6A3GtMUcQHLKqqn2hoTOHonJjIjLY0Hln5ERJSfSUNTWRGZxdvb0zln8FiuGTWJzC3lZG6pIDO7gq07Klm3uYJ1myvY+/0IIiNc6dhuKdEM2T2O+e/k0lNL6dU1kq7JEU1KRRLqy1o2p2JCJZgJ7pOAeFV9WUQGAa/iFuPtAM5V1XdD3ssg2AR3+7UxH/71FZyYChMH1H98oMZ8eBaX+snKriAru4J1W8rZuK2Czdsr97nzqiYRpdhXSFJn5Ycje9OzayQ9u0bSIyWCbsmRdEuKIDqq9giXXVhIv/vv3+eyFsDNixczJz2dTTfc0OhRwQUvvVRr8Al2Qj+wrzZSaZ9CNsGtqosCHmcCo0QkBcjT9rpow7Qq63fBNztd3Ytxvfad7K5PQnQ086ZPZ/aUKWTm5TGoS5c6P+w6xfpIGxRD2qB9i4KXlPnZ8n0lm7+vYN3WMl76cj3f71RiquKJ8Xcirqoz5Xnw+sdFtZ63c7yPbskRdEuOoHuyCyJdEn1sLc0lsaw7Y5NTKS71ExcjiEjQl7XCMafSVi6TWTALvWAuQ+0RUCkvNzTdMaZ+P+gPb2e5gPFcBlw+tvHn6J2YGPSHSFyMjyH9ohnSL5rjiefy01LILiwkMy+PPvHJ+Mrj2JZbybadlWzfWcnWnVXk5FWyY1cVO/OrKCjyU1Dk31M9cC/hME5l7r9hLpuJiRK6dPZR4K/i0KLJvLEommXJeSTF++hcvSVEkNjJR3yceH/6iPBWsodjTqW1XyZrK8GsLQoqWIjIL4Bf4ia3EZFs4D7gARtdmHATcQHi12/D4kxXJGl4t5btU83gM6BX7Xmr/H5lV6GfHflV5OzaG0DyCqrIK/Tz+focCoqUeOIpq/CxbWcVEEcKffn4q3KgvN6+xMcK8Z18xMZEMjZ3Kg8/U87n/XYQH+sjLlb4cEsB/YpHsWltJ97NKSI22kdcjBAX6yM2WoiLEWJjfMRECVGR7tZkCM9IJdTBp7UHs2ptcSRllfJMm9Q/Cc4YDi9/A498CX+Z5NZgtHY+n5CSFEFKUgTD+u//YbO7PMl9OK3MQPwRxGgnpg4cybVjJ1Be5nOjkt1+8ov8FBRVUbDbz+4Sbyv2U1Sq3ubmVZLpSf52WLK9OKCVbgyhG4++WEzgxH2t/RWIiRZiogR8VYwrOpONH3XmuhXb9wSTwrKhDN8Vwf3z8uiXXEF0lBAVuXeLjIDoSCGyel8EREQIBeUlvPllHlelTWJm6uFERgiXD+9PUX4MTyz/ii/G7qJX5wQifLgtQvCJq7zo8+E2ESK8x1t37271wWxP8FmZgSr4xMePRo7iwamnEBcZjd+v+P3gV/fFosoPfj9U+ZWqKvCr+7PKe62wtIJ7PvyY99ZvoMRXRGlUQdhuYAhmgjsXuKpGpTxE5BzgEVXtGsL+Bc0muNu/iiq45W3YUghnj4AZaS3do9Cpvqx1oDmV2lT5lZJSZXeJn8JiPzsLy7nn/U9I37ANnz+SKI1mdLfenNBvCBUVQkmZn9JypaRM3eMy97jU21/lr7/N1kJQqlAifUJkhA+fuFFolfopqiinc0wMUREuZ5iIWwMgAux5vO+NB1V+P9uLikiIjqZzTAyqigKFZeUUV1TQLa4TPhEUd3ee23Tf54D63Yd89Ye+ay30xh1WhQxY2+QbGEK5ghusUp5pBaIi4KrD4Pb34L/fwvi+bsTRHgQ7pxLhExI6CQmdfPTsCkOJ5uVRPww6+FRWKWXlSlmF+/O61xbx5tr1XDDqUA7t0Ysvt2znv9+sYULvAVw+ZhzlFUp5hVJRqVRUKRUVSkWVy+VVWVm93523qKySJZnrGdA5mYFJKVRUKZVVypb8QnYUlTK4SwqiQpUft1VVf9NWqtR94w78Jq4IPgS/H8r9gV+ChShiKCmFEhoX/WLoREUp7CwNvAMukmgiKShSXDhojL2Bwo2MXDCr1CoSoqO8UVP1aGnvqKn6cfWoKjJCqNIqvtyWTd+kRIZ1S+GYYXGccVxq2OrTW6U806aN6AZTBsGbmfDIF/CHE9x/QLOvYINPZIQQGSfEx7nn/z5/MlctXMjcjP9v79yjparuO/75StCIl/AQFF8IPopoY4wiEaPGaFWkGqtJUNSK2pqyDCu6NE1tNVaTrCbR0KyYRiOiCSpGVmKt1KpIEklcgvgCUXzgVa6CAhofiGCiwK9/7D3cw3Bm5szMmbkD9/dZa9Y5Zz9/57fP2XvOfv1+z4bFsQ//kAOYfPLhNXV7nHnXHKYvmsWlB45iTLSCePPcuZx+aHX/jM2MDQZn33U3v37uOS7+3GF8Yc8hzF7SwU8fe5zT9hvOLcbLcAAADs9JREFUjSedHMOy8d8/G883r/RXrlnDiMmTmXDICC4/6kggfH18749/4IYnnmT+P32NQb3bgM5nTgpdYyS+XArdZY8sfY1jbp3K/Wedyeh9O9cvz2xvZ/S0aTx83nkcMXiPzPf88KuvMuWX9/EfZ53FCfvstNG9UfbpMzUWaZby4nqLzSzl5SaZ42TkzL+GJ5dD+7twfzv87b5dLdHWSzVTj7Mw+eRQgU+aO5dr58zZZAC5GiTRQ3DTKSehbYwfPzaHSfMeCel9+gAmnzyGtioHtQb07cOpB+7FdQsepuf26zaa9L1u/lxOP/AAhu9a3WfssIH92WYbeKhjySaNRa2LOPfu33+jffoT9qk/vUpkNav6UMb0zMyOqVoIaTTwE6AHMMXMflDkvx3hi+YQ4G3gdDPrKJemj1l0L55cDtfMge16wLXHwc7lNo9yWo5au8kanV7es6EKiyQvHTVqM3vytYwx5J0etLDxI0k9gMXAcYT9pR4HxpnZc4kwFwIHmtkESWcAp5rZ6eXS9cai+/GTeTBnGRy5B0wc2dXSOFsTrdr4NGJqb97Gj3YmjFvsT+jyWwRcb2Zv1pDWKOAqMzshXv8rgJl9PxFmZgwzNxpeWgEMLLemwxuL7seqv8CMF+Erw8M+Uoveglkvw/4D4fi9Q5g/rYXb06ZnlCEt/o694O8T+ytf91gYcM1KqfgTR8In4jSR6YuC/Y5qSIs/9oCway/AQx3w9Irq0kyLf/QQOGhQcCvouRrS4nencjpt+GpeWxUan4df6113OQ3us5addvhTLl9muc2GkvR54AFgJZ3rLM4GLpF0gpnNLRk5nd2ApYnrZYQxkNQwZrZO0ipgR8J+VEnZvkZYA8LgwYOrFMPZ0umz3aYv9ltrYO7rYf1FoRJZ+3Fwq4a0+Hv03jSvectgXRX/u0rF/3oijWffhMVV7o2QFn9MYgyn473q7z8t/rABwZANdOq5GtLid6tyOrQ3g/v03iR+PeXUb/tenDa8sXVeLbOhfgT8CphgZhsAJG0D/ByYBBxeZXppc1eKizNLGMxsMjAZwpdFlXI4WxnDB8JFI2FgYvxiQK/gVg1p8bcvWqA9cWTnTrhZKBW/R2Ly+dj9YXXlBdubkBZ/UFun29F7wrAqV0KlxR/at9OtoOdqSIvv5dTpVm057dJWOUy91LIo70PgIDN7sch9P2C+mW1fZXreDeU4jtMilOqGqmUR3SqC/e1ihgLv1ZDe48C+koZK2hY4A5hRFGYGMD6efwX4ve9B5TiO0zxq6Ya6E7g5mlGdQ+gOOgL4AaF7qiriGMREYCZh6uwtZrZI0neAJ8xsBnAzcJukduAdQoPiOI7jNIlaGotvEcYQbknE/xi4AbisFiHM7D7gviK3KxPnfwa+WkvajuM4Tv3UYvzoI+CiOLawN6HhaDez8ttXOo7jOFssVY1ZSOopaZ6kYWa21syeMbOF3lA4juNs3VTVWJjZx4SBbB9cdhzH6UbUMhtqKnBB3oI4juM4rUst6yyuB84ClgBPAptYpTezb+QmXR1Iegt4NYekBlC0UrwFaEWZoDXlcpmy04pyuUzZyUuuPc1sYLFjLY1FuR1oa9p1tpWR9ETaApWupBVlgtaUy2XKTivK5TJlp9Fy1TIb6ouNEMRxHMdpXTKPWUjqJelnkl6X9KakOyQNaKRwjuM4TmtQzQD31cC5wP8RVnEfR1iIt7UzuasFSKEVZYLWlMtlyk4ryuUyZaehcmUes5D0MnC5md0Zr0cCjwCfNLP1ZSM7juM4WzTVNBYfAUPN7PWE24fAX5nZ0tIxHcdxnC2darqhegDFu7avo7b9pVoOSaMlvSipXdJme1xJ2k7S9Og/T9KQBsuzh6SHJD0vaZGki1LCHC1plaQF8XdlWloNkK1D0jMxz832gVfguqirhZIObrA8wxI6WCDpfUkXF4VpuK4k3RLH855NuPWXNEvSS/HYr0Tc8THMS5LGp4XJWa5rJb0Qy+duSX1LxC1b1jnLdFUcEy2U0ZgSccu+qznLND0hT4ekBSXiNkpPqfVAlzxXZpbpB2wg7Aw7I/H7GPht0i1req30IzSELwN7AdsCTwP7F4W5EPh5PD8DmN5gmXYBDo7nvQl2yotlOhq4twv01QEMKOM/BrifsG/YYcC8JpflCsJc8abqCjgKOBh4NuF2DXBZPL8M+GFKvP7AK/HYL573a7BcxwOfiOc/TJMrS1nnLNNVwDczlG/ZdzVPmYr8JwFXNllPqfVAVzxX1XxZTAXeAN5O/G4nmDtNum2JjCRshviKhY0S7wROKQpzCkEHAL8BjpWUZsEvF8xsuZk9Fc9XA88TzMtuCZwC3GqBR4G+knZpUt7HAi+bWR4LMqvCzP5I2EI/SfK5mQr8XUrUE4BZZvaOmb0LzAJGN1IuM3vQzNbFy0eB3fPKr1aZMpLlXc1dpviuj6UGMwx1ylSqHmj6c5W5C8nMzqs1ky2A3OyAN4LY5fVZYF6K9yhJTxMa8m+a2aJGy0PYG+xBSQbcaMGcbZI0fe4GLG+CbGdQ+oXuCl3tbGbLIbz4knZKCVNKX83ifGB6Cb9KZZ03EyWdAzwBXBoruSRZ3tVGcCSw0sxeKuHfcD0V1QNNf65q2RtqayQ3O+B5I6kNuAu42MzeL/J+itDd8hngp8D/NFqeyOfN7GDgRODrko4q8u8qXW0LfAn4dYp3V+kqC12iLwBJlxPGHqeVCFKprPPkBoLZg4MIfywmpYTpKl2No/xXRUP1VKEeKBktxa1mXXljEVgG7JG43p3w7zM1jIId8D7U9hmdGUk9CQ/INDP772J/M3vfzD6I5/cBPdWEhZJm9kY8vgncTegaSJJFn43gROApM1tZ7NFVugJWFrrg4vHNlDBdoq844HkScJbFTu5iMpR1bpjZSjNbb2YbgJtK5NV0XcX3/TRKf301VE8l6oGmP1feWARazg547CO9GXjezP6zRJhBhXEThXUv29DgcSNJO0jqXTgnDJQ+WxRsBnBOnBV1GLCq8MncYEr+++sKXUWSz8144J6UMDOB4yX1i7Najo9uDUPSaOBfgC9ZCXs0Gcs6T5mS41qnlsgry7uaN38DvGBmy9I8G6mnMvVA85+rvEfvt9QfYQbPYsJMi8uj23cILxPAJwndG+3AY8BeDZbnCMIn40JgQfyNASYAE2KYicAiwoyQR4HDm6CnvWJ+T8e8C7pKyiXgZ1GXzwAjmiBXL0Ll3yfh1lRdERqq5YRZgsuAfyCMa/0OeCke+8ewI4Apibjnx2erHTivCXK1E/qzC89WYabfrsB95cq6gTLdFp+XhYTKcJdimeL1Zu9qo2SK7r8sPEeJsM3SU6l6oOnPVdW7zjqO4zjdD++GchzHcSrijYXjOI5TEW8sHMdxnIp4Y+E4juNUxBsLx3EcpyLeWDiO4zgV8cbCcZytEkmDJX2Q+N3W1TJtyfg6C8dxHKci/mXhVEU0ujJTwZDQkw1I/35J38o73VZD0u6STA00otXVuuzq/J182Sqs3Dm1I2k28Fsz+17GKBOANmBH67SHkFveZnZiPWnWkmee4buKZutS0geJy+3i8S+JvNsaXZZZkPQg8ICV2F/NyY5/WTjVshdhU7O6GoqtibgraLciNgZtZtZGML4zrcitVfgsYXt6p068sXA2Eu0I/5uk38UBwWclHZ7w/1/CDpfjo//VknpJ+pGkJZLekfSApH0Scdqi/yuSVivYET5C0n8RDMp8O6b1Ygw/W9IV8XyipPlFMg6VtL7QfVMp/3ruuYyMHZKuVLCNvAb4sqQdJd0qabmkFZKmSuqfyGeQpBmx+24xKRbLYrpnJ66HxK6q3evVZbyuJGPZ8q+WlPw7JF0R9faBgs3qAyWNU7CnvUrSFIUtwesq23iva4EBwL0xv3+s9V4cbyyczTkf+AbBXscsOk03YmYnE4zkTI3/IP8dmALsR7C1PYhgxevexL/tmwmWzI4FPkUw/7jCzCYCDwPfjWkNS5FlGjBc0kEJt3OB2WbWEa8r5V/zPVeQ8QLgEkKX3D1R1n4E+8jDCZVUcvbNNGA9MJhg6/ncKuQrUI8uCzKUkxHKlH9OjCfYs+9H2KX1buCLwGeATxOMV42NYWsuWzN7G/gy0JH44pmS8710K7yxcIq50cwWmdl6wsu6j6Q+aQEVjAeNAy60YLjmI+BqgpH5zymYehxL2N55iQVeMrP2LIJYMKl5D3BezE+EyuaWLPk34p4T3GRm8y1MJ+xHsHd8iZm9G+W+BBgjaRdJuwHHEEy5rjKzFVHOzNSrS0m7lpMxEbQWXVTDZDN73sw+Bu4gdGtebmZrzOw1YDZwaE5lOwLIfRJGd8UHuJ1ikkaK1sRjb2BVStih8bgw1OMb6Umw0PVRvF5chzy/AG6X9M+Erpa+QMFaWKX8s1LNPRfoSJwX8lqScHs54VcQ7tWEfzJsFobEY626rCRjQQe16KIakumvBdab2VtFbr3Jp2xHEGyXODngjYVTD4XKb9+iFx7Y+G8YYF/guZT4GzLk8SDwZ4L5z1OBO83swyz550QpGZPuS+NxCMHIDIR/zAW/HvF8Tzor6EJlmOQDYIfE9a6J8454rFWXlWRsNfIo20MIRricHPBuKKdmLNgbvgO4Pna1IKmvpFMltUX/30T/IQrskxikXAGUHbC0YI/5VkI/+mnELqgs+ed0m1lkfIPQqE2K+fcDJgH3m9lyC+Y4ZwPXSPqUpJ2Bb6ck9QQwLg5kD0yGqVeXlWTMoohmkqVsJV0lqaNMMgPp/Kpz6sQbC6deLgBeBGZLWk0wi/lVgilICAOmC4A/AKsJYxCDot+PgRGS3pO0qEwevwC+ACwxs8eqzL9essp4NuH+Xoi/94BzEv5nEtYjLCUMRt+aksYVhEHw5YTG5c4i/3p1WUnGVqNS2Q4m6KkU3wd+FWdCHdNIQbsDvt2H4zhbJHEK8rFm1ordaFsd3lg4juM4FfFuKMdxHKci3lg4juM4FfHGwnEcx6mINxaO4zhORbyxcBzHcSrijYXjOI5TEW8sHMdxnIp4Y+E4juNUxBsLx3EcpyL/D0IH1dEwHnXZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.ticker as plticker\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(t_range, instant_rate, label = \"Instantenous Rate\", color = \"dodgerblue\", alpha = 0.75, linewidth = 2, linestyle = \"-.\")\n",
    "ax.plot(ode_time, 1 - ode_sol, label = \"ODE Approximation\", color = \"royalblue\", linewidth = 2)\n",
    "ax.scatter(list(final_size.keys()), list(final_size.values()), facecolor = \"white\", color = \"teal\", label = \"Simulation\", linewidth = 1.5)\n",
    "ax.legend(loc='lower left', bbox_to_anchor= (0.0, 1.01), ncol=2, \n",
    "            borderaxespad=0, frameon=False, fontsize = 12)\n",
    "ax.set_ylabel(\"Probability of Major Epidemic, $p(\\\\bar{t})$\", fontsize = 14)\n",
    "ax.set_xlabel(\"Infective Introduction Time, $\\\\bar{t}$\", fontsize = 13)\n",
    "\n",
    "\n",
    "#plt.savefig(\"imgs/exp_prob_epidemic_final.pdf\", bbox_inches = \"tight\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}