ExploringInfectionDynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this notebook is to explore the dynamics of the hoof and mouth disease based on the Keeling algorithm and the state space defined below. Of particular interest is the number of infected farms throughout the course of the outbreak. Insight into this particular disease transition will help determine when data for the risk model should be collected, and when an intervention should be implemented. NOTE: For now I am assuming an outbreak ends when all farms become infected (this usually occurs within 2000 timesteps, however for purposes of computational time I have truncated outbreak time to 1000 time steps). The code in the below cell block generates the location and size of farms. This state will be used to simulate 100 outbreaks and the number of  infected farms per time step will be averaged over the 100 simulations\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from random import randint\n",
    "Size = 20\n",
    "N = 75 \n",
    "np.random.seed(53080)\n",
    "x = Size*np.random.rand(N)\n",
    "np.random.seed(23003)\n",
    "y = Size*np.random.rand(N)  \n",
    "#np.random.seed(10)\n",
    "#Cows = np.array([randint(25,250) for p in range(N)])\n",
    "Cows = np.array([randint(25,51) for p in range(15)]+[randint(51,76) for p in range(30)]+[randint(76,95) for p in range(20)]+[randint(95,250) for p in range(10)])\n",
    "#np.random.seed(11)\n",
    "#Sheep = np.array([randint(25,250) for p in range(N)])\n",
    "Sheep = np.array([randint(25,51) for p in range(15)]+[randint(51,76) for p in range(30)]+[randint(76,95) for p in range(20)]+[randint(95,250) for p in range(10)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x106a23be0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEKCAYAAAAhPD1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4U9X/wPH3SZsmuUkHe6+yFdkgoEBFQFBEhgNwoIii\nqPh1o6CC/lzgF0G/gIoggooMURAQEKHIUJkCArIR2ZuOdCaf3x8JWqEjbZOmLef1PPd50txz7/nc\n2zaf3HPPPUeJCJqmaZqmBY8p2AFomqZp2pVOJ2NN0zRNCzKdjDVN0zQtyHQy1jRN07Qg08lY0zRN\n04JMJ2NN0zRNCzKdjLUiSyl1t1JqSbDjyIpSqrpSyq2UKpT/Z0qpF5VSk4JQb0+l1F9KqXilVKOC\nrl/TCiOlnzPWCjOl1PXAKOAqwAXsBP4jIhuCGpgPlFLVgf1AqIi4M1l/EHhQRH4sgFhigOkiUiXQ\ndfkQyz48v8PvsljvBpzAxQ+nNBEpWVDxaVowhAY7AE3LilIqAlgADAJmARagLZASzLj8SPgn4VwR\nlFIKqArsyKFoQxHZn8c6QkUkPS/balqwFMrmM03zqgOIiMwUj2QR+UFEtgEope5XSq3yvn7e2+x5\ncUlTSn3qXReplJqslDqqlDqslHo9q6ZjpVRLpdTPSqlz3vIfKKXMGda7lVKDlFK7vWX+l2GdSSn1\nrlLqlPfq75a8HLRSyqKUGquUOuJd3lNKhWVYf5tS6jel1AWl1F6l1E3e9x9QSu1QSsUppfYppR72\nvm8Hvgcqes9NnFKqglJqhFJqeob9dldKbfce1wqlVL0M6w4qpZ5RSm1RSp1XSn2llLJkEb9SSg33\nbnNCKfWZUirCWz4eCAG2KKX25OKc1FRKLVdKnfae38+VUpGXxPe8UmorEO8t7/b+jRxSSp1RSj2i\nlGqhlNrqPcYPMmxfSym10ntsp5RSX/kam6b5hYjoRS+FcgHCgdPAVKALUOKS9fcDqzLZrjJwBLjJ\n+/M3wETABpQBfgUezqLOpkBLPF9Uq+G5gnsyw3o3MB+IAKoAJzPU8wieZvRKQAlgBZ6mdVMWdR0A\nOmTy/mvAWqC0d1kDvOZd1xI4D9zo/bkiUNf7+maghvd1OyARaOL9uT3w1yX1vIqn6Ro8X3wSgBvx\nJMvngD14mtgvxvoLUN57bDuAQVkc1wDvttUBO/A1MO2Scxidze/dDdS85L2a3tjM3nOyEngvw/qD\nwCbvubd463YDE4AwoBOeFpVvvNtXBE4Abb3bzwBe9L4OA9oE++9fL1fWoq+MtUJLROKB6/E05U4C\nTiql5imlyma1jVLKBswDxorIEqVUOaAr8JSIJInIKWAs0CeLOjeJyDoRcYvIn8DHeBJZRm+LSJyI\n/IUn4V7shHQnngRxRETOAW8CKg+H3g9P8j0tIqeBkcC93nUPApPFe59ZRI6KyC7v60UicsD7+idg\nKZ5mfbKII+N7dwELRORHEXEB7+L58tImQ5n3ReS499i+AxpnEf/dwH9F5KCIJAIvAn1y2ZFtk/fq\n9ZxSaqyI7PPGluY9J+/x79+LeOM7IiIZb2O8LiKpIvIDnqvyL73n9SiwCmjiLZcKVFdKVfKWX5uL\nWDUt33Qy1go1EflDRB4QT8ejBniuaMZms8lkYKeIjPb+XA3P1dSxix/uwId4rpAvo5Sqo5RaoJQ6\nppS6ALwBlLqk2PEMr52Aw/u6AvBXhnWHcj7CTFUE/rxkPxW9rysD+7KIvatS6hdvk+w5PFfKl8ae\nXZ1/xysigudYKmUok/G4k/jnuC9VIZP4Q4FyPsYCniv6Et7lP0qpct6m8cPe38t0Lj+2vzLZz4lL\nYr7053Dv6+fxfDlZp5T6XSn1QC5i1bR808lYKzK8V4Cf4UnKl1FKDQVq4bl6vOgvPM2TpTJ8uEeK\nyDVZVDMRTxNsLRGJBIbh+//JMTydky6qmlXBHBzF08yacT9HvK//wnOM/+K9H/s1np7nZUWkBLCI\nf65+c+oodgTPF5eL+1N4muGPZFE+u/1lFn86/06EufUmnib/Bt7fy71c/nvJS2c4ARCREyLysIhU\nwtNhcIJSKjof8WparuhkrBVaSqm6SqmnlVKVvD9XAfoCP2dStivwBNArYzOliBzD01w7RikV7u1k\nVVMp1S6Lah14mjOd3g5Mj+YUJv8kvFnAEKVUJaVUCWCoD4cZppSyZlhC8dy/HK6UKq2UKg28Anzu\nLT8ZeEAp1cF7LJWUUnXx3OcMw3OP3e09H50z1HMCKKU8PdQzMxu4xbtfM/AMkIzn3nVWx52VGcBT\nyvOctQNPIv1KMnm8KxcceO6Bx3n/Hp7Lx74yUgBKqTuUUpW9753Hk6TzE6+m5YpOxlphFg9cC/yq\nlErAk4S34kkU8O9Hg+7E0zFnp/qnR/UE77r78CSqHcBZPImnfBZ1Povnnm0cnvvFX/HvK65Lr74y\nxjAJWAJsATbguVLN6WptEZ6m7ovLK8D/ebff6l02eN9DRNYDD+C5Z3oeiAWqeu+vD8HzheAsni8t\n8/4OUuQPPElyv1LqrFKqQsbYva0O9wAfAKfw9AS/VbJ+RCi7x7Km4GlG/gnPc9ZOPF+UMm6bnczW\nj8TTue4CnvvVvpxbX66UL5ZpDvyilIrHc96GiMhBH7bXNL/wadAP7xXJNKAsnj/ej0XkfaXUCGAg\nnn9e8PRGXBygWDVN0zStWPI1GZcHyovIb95mp41ADzxXI/EiMiawYWqapmla8eXTCFwichxvT0oR\nSVBKXXyWEvL26IamaZqmaV65vmesPOPtNsEzAADAE95ReSYrpaL8GJumaZqmXRFylYy9TdRz8IxI\nlIDnMZAaeB7+Pwb81+8RapqmaVox5/OsTd7HHRYA34vIZYMueK+Yv7v0+U2l1BU1EL6maZqmiUiu\nbuH6dGXsHQBgMrAjYyL2Ph5xUU9gWxZB6SWH5dVXXw16DEVl0edKnyd9rvR5KsxLXvg6heJ1eJ5B\n3KqU2ux97yWgr1KqMZ7HnQ7gGblG0zRN07Rc8LU39Woyv4r+3r/haJqmadqVR4/AVUjExMQEO4Qi\nQ58r3+jz5Dt9rnyjz1Pg+NyBK88VKCWBrkPTNE3TCgulFBKIDlyapmmapgWOTsaapmmaFmQ6GWua\npmlakOlkrGmapmlBppOxpmmapgWZTsaapmmaFmQ6GWuapmlakOlkrGmapmlBppOxpmmapgWZTsaa\npmmaFmQ6GWuapmlakOlkrGmapmlBppOxpmmapgWZTsaapmmaFmQ6GWuapmlakOlkrGmapmlBppOx\npmmapgWZTsaapmmaFmQ6GWuapmlakOlkrGmapmlBppOxpmmapgWZTsaapmmaFmQ6GWuapmlakOlk\nrGlakeRyuVi0aBExMbcSEVEOs9mKYURRo0ZDRo/+L2fPng12iJrmMyUiga1AKQl0HZqmXVmmTfuc\n5557haSkUsTHPwbcBEQBScAObLaPEfmO22+/g4kTx+BwOIIbsHZFUUohIipX2+hkrGlaUfLSSyMY\nN+5znM7PgVbZlDyF1foC1aptYdWqxZQpU6agQtSucAFLxkqpKsA0oCwgwMci8r5SqiQwE6gGHATu\nFJHzl2yrk7GmaX4xZsz7vPzyhzidsXg+jnIimM1DqVfvJ379dTk2my1gsZ09e5bp0z/nt9/+ID3d\nRd261ejf/16qVKkSsDq1wimQybg8UF5EflNKOYCNQA/gAeC0iIxSSr0AlBCRoZdsq5Oxpmn5duTI\nEWrVuobk5E1A9VxsKdhsvRg+vBUvvfSC3+NKSEhg8OBnmD17FibTLTidrQETFsvvwAzatYvh008/\noFKlSn6vWyuc8pKMferAJSLHReQ37+sEYCdQCegOfOYt9hmeBK1pWiGVnJzMtGnTaNCgDWFhBmaz\nlVq1mjBhwkQSEhKCHV62Jk6chEgfcpeIARRJScMYO/ZDXC6XX2OKi4ujZcsbmD07jeTk3d6m88eA\nR0lJGU9KyiFWrGhE48ZtOHjwoF/r1oqXXN8zVkpVB1YCDYBDIlLC+74Czl78OUN5fWWsaYXA4cOH\nadeuK6dOVSAh4XGgA57v42ux2ydgt29m5crvqVevXpAjvVxaWhply1bj/PklwDV52kd4eEtmzHiV\nW265xW9x3XprH374IYqUlIlA1hdCJtP71Kgxmd27N2My6YdYiruAXRlnqMABfA08KSLxGdd5M67O\nuppWCCUmJnL99Tdx6FBfEhKW4GnUcgAG0JHExLmcOvUy11/fmZMnTwY32EysX78el6sceU3EAPHx\n9/LFF9/4LaaDBw+ybNkyUlLGkF0iBnC7n+DECRPLli3zW/1a8RLqa0GllBlPIp4uIt963z6hlCov\nIseVUhWATP+LR4wY8ffrmJgYYmJi8hywpmm5N23adE6dqonL9VKWZUQGkJCwjnHjxvPGGyMLMLqc\nnT59GqUq53MvlTh+fLlf4gEYP/5j3O578XyhyYkiIWEwo0ZNoHPnzn6LQSscYmNjiY2Nzdc+fO3A\npfDcEz4jIk9leH+U9713lFJDgSjdgUvTCp/o6EYcODAGuDGHktuJiurEyZN/YjabCyI0n3z77bf0\n7z+FuLj5+dkL7dt/SmzsPL/E1Lp1V3755XHA12bvg5Qq1Z7Tp//0S/1a4RXIZurrgHuAG5RSm71L\nF+BtoJNSajeeG1Bv5ypiTdMCLiUlhT//3IHnXzQnV5OWFsrRo0cDHVaulCpVCjiez72coGzZkv4I\nB/CcV7DmYgsLaWkpfqtfK158aqYWkdVknbg7+i8cTdP8LT09HZMpFLfbty/qSoWSnp4e4Khyp0WL\nFrjdB4G9QK087SM8fAZ33vm432KqXLk8mzfvI+fWhov2U6pUeb/VrxUvulufphVzhmHgcJQAtvtQ\n+gTp6ecoX75wJQ2r1crAgQ8QFvZhHvewnZCQPdx2221+i+nhh/sRHj7Z5/I22yQefrif3+rXihed\njDWtmFNKMXjwQ1gsE3IsGxLyCXfccQd2u70AIsudJ54YhMk0lSz6iWbLYhnF4MED/XofvGvXrlit\nJ4AffSi9H5FvGThwgN/q14oXnYw17Qrw2GODCAubA/yQTalNWCzjeOGFJwsqrFyJjo7myScHYxi3\nAvE5lr8oJOS/VKiwkeeff9qv8YSEhPDFF5Ow2foB67IpeRDD6MI777xB6dKl/RqDVnzoZKxpV4CK\nFSuycOEc7Pa7CQl5jX93hjqHUmMwjC58/vnHXH311cEKM0dvvTWSO+9simG0A/blUDoJs3ko5cp9\nSGzsQiIjI/0eT6dOnZg58xMM42as1oHAJv4ZbmE/ZvNzWK3Nee21xxgy5DG/168VHzoZa9oVom3b\ntmzatJp+/Y5gtdYnIqIJERFNsVqj6dFjE6tXL6Fnz8I9oq1SiilTJvDyy/dgt1+Lw9ENWARc7HAm\nwB7Cwp7Faq1Ku3Y7+e23tVSrVi1gMd16663s37+dYcNqUrp0L0JCrISG2rDbW/Loo7Bt2y8880zh\nbG3QCg89haKmFTMnTpxg5syZHDx4mPPnEyhVKoLataO56667/r46jIuLY+/evbjdbqKjoylZ0n+P\n/BSUpKQkvvrqK955ZwJ79mwmJMSGy5WMw1GSBx7oz5AhjxAdHV3gcSUnJ+NyuTAMA88QDdqVRs9n\nrGlXKBFh9erVvPvuBJYuXYxSPUhKqgfYgXgMYxNu9zLuvPNOnn56MI0aNQJg48aNjB37ERs2bMVq\ntdKnzy08+OADRe7eptvtJj4+HpvNRlhYWLDD0a5wOhlrOTpx4gTbtm0jNTWVKlWq0KBBA/3tvYhL\nS0tjwIDBfPPNcpzOIYj0B6IyKXmMkJBPCAubwNNPD+LcuTimTp1NcvJg3O52QAI22wxCQhaxYMFs\n2rdvX8BHouUkLS2N5cuXc+zYMcxmM40aNaJBgwbBDku7hE7GWpY2btzIiBGjWbZsCRZLEyCM9PQ/\nqFixFC+++Dj9+/fXs8kUQS6Xi5tvvp3Vq1NwOmfhmfwhJ8cJC7sJl+sMLtdW4NIm6h9xOPqyY8dG\nqlSp4v+gtVyLj4/n7bffZcKESbjd1XG5aqNUKm73KmrXrsGrrz5Nz549gx2m5pWXZIyIBHTxVKEF\n05w5X4thlBGlxgqcExDv4hL4XgyjufTufY+kp6cHO1QtlwYNelIMo6NASobfqy/LeYGrBMZmuj4s\n7Al5/vlhwT48TUROnToldeo0Eau1r8Dvl/yuUgW+FsOoJc89NzzYoWpe3ryXq1ypr4yLuY0bN9K2\nbVeSkpYATbIolYRh3MIjj7Tiv/99syDD0/Lhzz//pF69piQn7yPzZumc7AWuBQ7hubec0QYqV76P\nv/7akd8wtXxwuVw0a9aOHTuuJy3tbbKeqvE0htGOUaOG8NhjjxRkiFomAj6fsVb0vPLKOyQnv0zW\niRjAhtP5JRMmTOTChQsFFZqWT/9M4ZeXRAyeMZ6vA2Zksq40SUmJeY5N848lS5awb19yDokYoDRO\n5wyGD3+dtLS0ggpP8yOdjIuxY8eO8eOPP3g79OSkPCZTFz77bFrA49LyLyUlhQ8//ITU1EfzuafB\nwHj+GajionVER9fO5761/Bo1agIJCY+RfSK+qBEuV03mz8/PNJNasOhkXIxt2bIFq7U5EOFTeafz\nJmJj1wc2KM0v1q5di1I1gLr53FNn4BiQcY7ddOz293j22YfzuW8tP0SENWuWAbf7vE18/B3Mm7c0\ncEFpAaOTcTGWlpaGiCUXW4SRkpIasHg0/zl9+jQilf2wJxNQCTjt/fkkVus9NG4cTq9evfywfy2v\nXC4XLlcaEJ6LraI4fz4hUCFpAaSTcTFWuXJlXK4/uLwJMnMhIX9Qs6Y/PuC1QEtJSUHEX4NbpGCz\nvUBExM1YrXXp27cUP/wwj9BQn6Y71wIkNDQUs9kKnMvFVqcoUyavfQi0YNLJuBhr3LgxZcoYwHIf\nSqcRFjaZhx7y5f6yFmxRUVGYTOf9si+7PY2RI7swdepDHD26nylTxmOz2fyyby1/OnXqhlJf+lha\ncDi+4I47bg1oTFpg6EebirmPPvqYp5+eitO5Asi6ydpkGkOTJvPZsCG2wGLT8s7zWFMzkpMPAUY+\n9nQIm60xp08fxjDysx8tEFauXEm3bo+QkLAVyGku5tWUK9efo0f36AF8gkw/2qRd5qGHBtKhQxUM\noztwKpMS6ZhMY4iKGsOcOVMLODotr6pVq8a117YCvsrXfkJDP+bee+/JMRE7nU6mTJlCw4bXU7Jk\nFcqVq8ntt9/HL7/8gv6yHTjt2rWjRYuaWCwPAa5sSh7AMO5m3Li3dCIuovSV8RUgPT2dJ598gSlT\nJqPUrSQl3QSEERKyA4vlE+rWjWbu3GlUr1492KHm25kzZ5g6dRrr1m0jJSWVqlXL079/X5o1axbs\n0Pzu+++/5847h5OQsAHfHn25VCo2W1U2blxB/fr1syy1ZcsWOna8leTkRiQkDAIaAsmYTN9hs02g\nY8eWzJw5FYslN50FNV8lJCRw0009+e03F07nUKAj/1xHnUepadhsb/P228N54onBQYxUu0iPTa1l\n68yZM0yZMpWfftpISkoqtWpVYdCg+/+ewScQTp48ybJly4iPj6dUqVJ07tyZiAjfHrXKjbi4OB59\n9Bnmzp2DUt1JSmoLhGEyHcBq/ZRq1coxadIYrrvuOr/XHSxut5saNa7mr7+eR+SBXG9vNr9Mixbr\nWLNmSZZlDhw4QJMmbbhwYRxwZyYlkrHZ+nDTTXa++eaLXMeg+SYtLY1PP53KO++M58SJC4SG1gRS\nSU7eRqdON/HSS0/SunXrYIepeelkrBUaBw4c4Omnh/P99wsJC+uAy1WKkJAjpKf/TJ8+fXn33f/z\n2xy658+fp2XLGzh0qDkpKW8Bl07/5wK+xjAeZ9asT7nlllv8Um9h8Mcff9CyZXvi4z8Cevi8ncn0\nP8qWfY/ffltLuXLlsix3770P8+WXFXC7R2azt2Ts9vqsWDGLFi1a+B68lmsiwvbt2zl69CgWi4V6\n9epl+/vTgkMnY61Q2LFjB9dd15G4uEdwu58ASmRYewSz+Q0qVFjBunWxfvkgiYm5hZ9/rkVq6liy\nb679FcPoxm+/raV27eIzutSGDRvo1Kk7cXGP43Y/TvaDvJzEbH6DsmW/56efFhMdHZ1lyfPnz1Oh\nQg2Sk/8Asv89mUyjuP32ncyc+WmejqEwERG2bt3K4cOHsVgsNGvWjBIlSuS8oaZ56VmbrjBut1v2\n798vO3bskMTExGCHIyIiqampUr58tCj1WbazBpnNw+Xaazvku74tW7aIYVTyzl6T82xFoaHDZNCg\nIX440sJl//790rXr7WKxRInF8ojAJgGngFsgQWCVGEY/sVojpW/fAXLq1Kkc9xkbGysREW18nAVq\nh5QvX7sAjjRw3G63TJ48RaKjG4ndHi2RkV0kMrK9WK1R0rfvANm7d2+wQ8wzt9st69atkzlz5sjc\nuXNl27ZtwQ6pWCMPszbpZFxETZs2XapVayA2WwUJD68rNlsJeeihx+Xs2bNBjWvWrFkSHt7ehw/v\nNDGMyrJp06Z81Xf//Y9ISMjIXEwd+KcYRkmJj4/30xEXLkeOHJHhw0dI2bLREhpqEaVMYjbbpEqV\nq+Tdd8fImTNnfN7XDz/8IJGRN/h4XvdK6dLVA3hkgeV2u+XuuweK3d5UYKl3etGLx3ZCQkJGSERE\nOdmwYUOwQ82VlJQUef/9D6RKlfricNSRiIgeEhHRXQyjilx11bUybdo0cblcwQ6z2NHJ+Arxyiuv\ni2HUE/jRe+UjAockLGyQREdfI+fOnQtabK1adRaY4dMHeEjIa/Lgg4/lq75KleoLbM1FMhaJiGgp\nq1ev9tMRF15ut1tSU1PzvP2ePXvEZivvY6vDXGnaNMaP0ResN98cLYbRytuKkPUxlihRMaj/X7lx\n4cIFadnyBjGMTgKxGT4rRCBdYJ7Y7c2lZ89++fo70S6Xl2SsH0grYvbs2cPo0e97B/HowD/3SKuQ\nmjqRw4ebM3LkW0GL78CB/UBTn8q6XE3ZsWNfvupzOuPxdSKMi5SKJCGh+I/fq5TCbM5poIis1apV\ni3r16gDzcizrcEzk6acH5rmuYEpLS2PUqDE4nZ9w+bzOGfUkJaVtkZjZzO12063bnWzZEo3T+T3Q\nnn/3pwgBupOY+BNLlpzh4YeHBCdQ7W86GRcx48dPIj19AFA+k7WK1NQXmTz506DNaRoSEkL2gxNk\n5Mr3+Mfh4VH8M8mBb0ROERkZma96CwuXy8XChQt58slneeqp5/nxxx8vtkj5xWuvPYdhPI9nZqfM\nKTUdu30vt9/u++xChcnChQtxuWoDV+dY1ukczHvvTQp8UPn0ww8/sHnzMVJSPsSTeLNiw+mczcyZ\nc9m9e3dBhadlQifjImbTpp2kpbXJpkRt0tMVZ86cKbCYMmrWrAlK/ehTWYtlGW3aNM5XfT17dsVs\nzs0oVDuBYzRt6tvVe2H2559/Urt2Y/r2HcH775dm7NhIevT4D9dc04oTJ074pY5u3brx4ouDMIzr\ngJlAxlm9DhMaOpSoqKEsX76gyA76sWvXLpKSWvpYuiWHDxf+pOWZB/lxwJcvu+Gkpw9g3LgPAx2W\nlh1f2rKBKcAJYFuG90YAh4HN3qVLFtsWRBP9FeOWW/oITM7mvlaSmM12uXDhQlDiW7lypTgc9S/p\nAJPZEidWa0k5dOhQvurbu3evWK2lBRJ9ul8cFva4PP/8MD8dbfCkp6dL9epXicn07iXH6JbQ0GHS\nqFEbcbvdfqtv4cKF0rz5DWKzlZXIyLYSEdFCbLYS8uCDj+X7dxhso0ePFrP5aR/7HCRKaKg12CFn\nKzk5WUJCLDnc/7502SklS1YJdujFBoHqwAW0BZpckoxfBZ72YdvAH/kVZNasWeJwtMvmn2qatGrV\nKWjxud1uadEiRsLCnrykw0jGJUVsth5y990D/VLnHXf0F5vtdoG0HD5w5klkZHk5cuSIX+oNpnnz\n5kl4eKssjtMldnsdWbVqld/r3b9/v6xYsULWrFkjcXFxft9/MCxdulQcjsbZ/L1mXL6T+vWvDXbI\n2Tpx4oT3C6qvidjz5dhsNoIderGRl2TsUzO1iKwi80k18zIgrpYPPXr0oFSp04SEjMpk7U5sthd4\n442hBR7XRUopFi/+mtq1f8Zmuw1YBX/Pp+wCFmC3x3D99cKUKeP9Uue0aR/SvHk8NtutwLZMSlxA\nqVGEhz/E0qXzqFixol/qDaZFi34kPr53FmtNOJ29WbbMt9sFuVGjRg1iYmJo06YN4eG5mfS+8Lrx\nxhtxOOKBn3Ms63CM5/nnH/F7DG63m127dvHll18yatQo3nnnHaZNm8b27dtxuXztg+FhtVpJT3eC\nj/OYeyQRFuaZNjMtLY1t27axZs0atm3bFrT+J1ccX7M2UJ3Lr4wPAluAyUBUFtsVwPeQK8uhQ4ek\nZs2GEh7eXGCswKdis90nVmuUTJ06LdjhiYhIYmKijB07TipVqiuGUVkiIhqKzVZW6tdvKVOnTpX0\n9HS/1peamiojRvyfREVVlPDwtqLUSwIjxGbrLxZLlHTrdpfs2rXLr3UG06OPPikwKssrHZPpBRkx\nYmSwwywyPv30MzGM2gLHsjmnY6VKlbridDr9Vu+FCxfk7bdHSbly0WK3V5Pw8NslNPRZCQ19ThyO\nPuJw1JISJSrJyy+P9GmgFhFP61SlSnUFfsrFlfFUufbaDjJ8+Ajv/1BdiYxsLeHhdSUqqqIMG/aq\nHD161G/HXdyRhytjn4fDVEpVB74TkWu8P5flnzn5XgcqiMiDmWwnr7766t8/x8TEEBMTk8uvDNql\nXC4XS5Ys4auv5pGYmEyrVg0ZMOB+SpUqFezQ/kVE2LdvH/Hx8ZQsWZJq1aoFtL60tDQWLFjA77//\nTnJyKhUqlKN3795UqFAhoPUWtIULF9K37wji49dxeQOVG7u9DsuWfU6rVq2CEV6RNHLkm4wa9RFO\n50vA3YDDu2YTVus4SpRYw5o1P1CjRg2/1Lds2TL69n2QxMQ2JCX9B2hJ5o2NW7FaP8BiWcDkyf+j\nd++sWkSr2HqIAAAgAElEQVT+MXbsOIYN+wWnc4ZPsRjGVcBJXK47SEkZDFyTYe3vWCwTsFjm8v33\nc2nTJrsOpFem2NhYYmNj//555MiRSKCGw+SSK+NcrAvkFxCtCHC73bJ69Wq57bZ+Eh5eVsxmQ0qU\nqCQPPPCoHpYvj9LT06VmzYYSEvJ/l9zrdInZ/Iw0a9bOrx24rhQrVqyQTp16isUSIeHhdcVuryKl\nSlWR1157I1ejl2UlNTVVZs2aJdWrXy1QUmBxLq5e14hhRMvw4Tm3eJw7d04iI8sJLMxxv0q9IkrZ\nBRblUPZ7sdvLyNatW/N9Hoo7AtWBSzJvpq6Q4fVTwJdZbBf4I9cKrQsXLki7dl3Fbq8tSr0n8JdA\nnMAeCQl5VWy2inLvvQ/pEYDy4PDhw1K3blNxOBoIvCpKDRe7vY40bdrW5ybNwsbtdktsbKzccUd/\nadSovdSq1UyaNLlB7r//Edm4cWOBxXH27FnZsWOH7Nu3z2+3VBYtWiRRURXEYqkvUElgfy47WYnA\ncTGMejJmzPs51rdmzRqx28sITBPPiFuX7itFTKbREhJSUmC8j/VPDGoH0aIiYMkYmAEcxfOQ4V/A\nAGAasBXPPeNvgXJZbFtAh6/5y7Zt22TAgMESEVFOQkMtEhlZQQYNGiI7d+7M1X6Sk5OlWbN2YrE8\nJFn3dI4Xw7hJ7rjjPn0llwdut1uWLVsmQ4e+JC+9NFxWr15dJM+jy+WS8eMnStWqV4nDUV+UGiuw\nXOBXgWUSEvJ/YhhV5OqrW8nMmTODHW6uzZnztdhs5QTmCpQR2JiHRHxx2Sc2Wymf/h83bdok9eo1\nF8OoLibTGwLfCsyV0NAXxWYrJ40atRaLpYxAko91J4vNVrZY9b8IhIBeGed10cm4aBk58k2x2cpJ\naOgI7zd3p8BeCQ0dLjZbGRk9+j2f9zVmzFgxjK6S8zPHiWK3XyWLFi0K4JFphVVSUpJ063aHGEZr\nbwLO6hGjNIFvxTDqyeDBTxWZCQ4OHjwohlHam4DvExiWj0TsWUymcXLddTf5HMO6detkwIDB0rZt\nN4mJ6S6PP/6M7NixQ4YNe0VCQ5/JVd1m87Py4osvB/CMFX06GWv58sEHE7wTUBzJ4h/xTzGMmjJl\nytQc9+VyuaRixToCq3z8J58sMTHdCuAotcIkPT1dunTpJTZbr1xcnZ0Vw2gjjz32dLDD98kzzwyV\nsLCnBE4KRAqczncyhiSxWsvI7t278xVbv34DBT7KZd0fS58+D/rp7BRPeUnGejhMDYCUlBReemkE\nTufXQFbP4VbF6ZzJs88OJz09Pdv9bdy4kbg4BVznYwR9WLv2J86dy+xxdq24GjVqDD/9dJqkpC8B\nq49blcDpXMDUqYuYM2dOIMPLt5SUFD7+eAqpqY/guZvXFfDHEw9WXK5+zJgxM197CQ3NzVjyF6Vj\nNudvTHntcjoZawB8/fXXiDQErsqhZDNSU6uwYMGCbEsdP34ck6kmvo8LY2CxlOPkyZM+lteKuvT0\ndN59932cznFAbse1LkFi4pu8/vrYQITmNwcOHAAigTrAenz/cpqztLQ2xMZuyNc+rrmmDjbbr7na\nxmRahcUiuN3ufNWt/ZtOxhoAP/+8kYSEzj6VTUi4iQ0bNmVbxmKxoFRKrmJwu1OK7GQDWu4tWLCA\ntLQqQF4nC7mVvXv/ZMuWLf4My68SExNR6uK0jLvwZWYo313Nrl1/5GsP/fvfi9v9LeDrxDJncbsX\nMGPGaho2bM3x48fzVb/2D52MNQDS0134NsMLQIi3fNYaNmxIcvJG4IKP+9xJaGgqlSpV8rH8lePc\nuXOsXbuWH3/8kZ9//pm4uLhgh+QXY8ZMIj7+0XzsIZSUlIf54IPCO6VhiRIlSE8/BQiQgu9N8b6w\nkZqauy+8lypTpgzdunXHbP4/H7d4A+hBYuLv7NrVhVatOuhbS36ik7EGwDXX1MUwfvGprMPxK/Xr\n18m2TPny5enYsTNKTfdpnxbLRB599CHMZrNP5Yu748ePM2LE/1GpUl3KlatK167/oXfvN+nSZQil\nS1egSpWreOutUZw+nbu5nAuT3bt3Adfmax8uV0u2bt3ln4ACoEaNGpQo4cAzRnsEcNaPez9LeHj+\n5+X+8MMxlCu3iNDQ1yHL8awF+D/gO2AMoEhPH8mxY20YMeLNfMegoXtTax7nzp0TqzVK4GgOPSkP\niGGUlMTExBz3+csvv3ifrfwjh32uEIejjBw+fLgAjrRwS09Pl9Gjx4jNVlKs1ocF1svlAzakCawV\nm+1+MYxSMn78xCL5bHF4eFkf/t5yWtZK3botg30o2Ro7dpwYxp0Czwr8nx96Uv/Tq7l37/v8EuPR\no0elfv3m3tmrPhY4J55HzM4LTBJoItBMLn/SYq84HKX9Ol53cYB+tEnLjyeffEEMo6N4ni3O7J8/\nXgzjehk2zPcJCKZMmSo2W3nxDHZw6cAfToGPxG4vI8uXLw/gkfnXoUOHZP369bJy5UrZvHmznD59\n2i/7TUhIkOuv7yx2e1uBPT5+IG8Xu72FdO7cQ5KSkvwSR0EpWzY6F8eZ1bJUmjW7MdiHkq3z589L\nmTLVBJ4QiPFbMjaMO2XixIl+i9PlcsnixYvFai0vYBVQAjaBW8UzbGfm4wVERHSQb7/91m9xFAc6\nGWv5kpaWJj179hO7vanAHIFU7z9cisAMsduvkX79Hsz1YAtLly6VBg1ai2FUEbP5SYFXxWIZJFZr\nabn++q6yadOmAB2R/yQlJcn06dOlQYPWYrWWkYiIJhIZeZ1ERFwjFkukdO16uyxfvjzPV6ipqanS\nqtWNYrX2z+RLS05Lithst0vHjt2LzEAYIiLNm3fw/p3lPSEpNVr69Sv8z7z+/vvvYjZHCpQQ2OGH\nZHxMrNYoOXfunN9jjYysIHBYMh9CM7MvBQPk448/9nscRZlOxlq+uVwumTFjhjRu3FbCwiLF4YiW\nsLAIad78BpkzZ06+mkM3btwo77zzjrz88isyZswY2bt3r9/ivnDhgmzcuFF27Njh9ybb998fL+Hh\nZcTh6Cye4QQvTZYXRKnx4nBcLVWq1JOVK1fmuo5XX31dDKOzzx+AmSVku72NvPdezmMWFxaff/65\nOBwd85GQXGK315K1a9cG+1By9PXXX4vNVkfgaoEOkvUoY74tVusAefjhIQGJtUyZGgK7fY7F4egr\n06YVjqlbCwudjDW/OnnypOzevbtQTzpw+vRpueeeh8RqjZKIiEZit1eXSpXqyuTJn+Z73263WwYP\nfkoMo76PVzNugflis5WRr77yffzkPXv2iM1WWuBQPq+WdolhlCoy996Tk5O9941z6lOQ1bJEatZs\nXCTulzdvfoPAV+JpZbpKPPdl8/p7nidlylSTuLi4gMR62239RKkxPsbin5HAihudjLUrytmzZ6V6\n9avFbH5C4HiGhPiTGEZ9GTHijXzt/6WXXhXDaC5wNpcfllvEZisnS5cu9amewYOfErP5xXwmYs9i\nsQyWYcNezddxF6Tnnx8uNlt3yX2LQKIYRnP55JPJwT6EHJ04cUIslihvIhaB7QLlBGbn4Xf8oxhG\n6YC2BqxevVrs9tqS85jyIjBNWrfuHLBYiiqdjLUrynPPvSQWy/1ZfEgcFau1pPz555952ndsbKwY\nRtUMST63y0/icJSW8+fPZ1tPenq62O0lBQ74JRnDVilVqkqejjkYkpOTpUWLGAkLG5SLhOwUm+1W\n6d37niJxVbxjxw4JD69zyTFsFs80ikMEEnw45hQxm1+W8PAyeboNkhtut1uuvbaDhIX9R7JvTt8l\nNlt5WbFiRUDjKYrykoz1c8ZakeR2u/noo8mkpLyQRYkKuN138/HHU/K0/7fe+gCncyhQLo8RtkXk\nBqZNy/456927d6NUFJ7pwv2hAYmJTo4ePeqn/QWWxWLhhx++pVGjXdhsvYCdOWyxHsO4kc6d7Xz5\n5WSU8nW41eAJDw8nLe0skHH4yMZ4ZqA9C9QCXgG28+9xot3AbkymtzCM2rRvv4Vdu7bQrl27gMar\nlGLRotlER6/Cau3njSujJGAqNlt7/ve/t4iJiQloPFcKnYy1oLpw4QK//vorP//8M2fO+DokHyQk\nJJCUlADUy7JMamoztm/fl+uYjhw5wsqVy4F7cr1tRomJjzF69ISLLUSZ2rp1K0rldTjIzCgsliaF\neojIS0VGRrJq1WL+85/GRETcQHj4jcBs4A/gMLAD+JTw8BaUKXMnr79+F3PnfkFYWFhQ4/ZVpUqV\nqFy5CvDDJWtKAtOBH4FzwG1ABCEhNYiIaEBYWElKl+5Mnz77+OmnufzwwzwqVKhQIDGXLFmS9etj\neeaZukRFdSIiog0Oxz2Eh/fGaq1Ku3azWbx4JgMG3F8g8VwRcnspndsF3UytZWLfvn3Sr9+DYrVG\nSkREc4mMbCkWS6T06nWP7NixI8ft09LSxGy2CZzJshlNqTfl4YefyHVsw4a9KhbLYD80GbvF4bgq\n22a8Tz75ROz2AX5qovYs4eG9ZeZM3zuQFSYpKSkyY8YMufbaTlKhQh2JiqooFSvWkw4dusuCBQsk\nPT092CHmyaRJk8RubyfZP7YWL4ZRSyZNmiTbtm0rNB0nU1NTZenSpTJt2jSZOXOm7N+/P9ghFXro\ne8ZaUbB582aJiCgnJtMrl9yTPS1KvS12exlZvXp1jvvp2fNuUWp0Fh9saXl+7OX6628RmOeXxGg2\nPyWjR4/Osq7p06eLw9HHr8k4IuJmmTdvXq6PWwuc1NRUadu2i9hsdwrEZ/J7OyWG0V7uvvvBInEf\nXMteXpKxbqbWClRKSgqdO99GXNz7uN0j+fc92VKIvEBi4nRuvrl3jhMivPbaUGy2d4DFl6xJxmJ5\nkGbNatOqVatcx3j27Dk8TYj5l5ZWgtOnsx5Iv379+phMv/ulrotcrt+56qqcpsLUCpLZbGbJkrl0\n62bDaq1GWNiTwAzgC2y2AVittXnwwZZ89tlHReI+uOZ/eoZorUDNnTuXpKRawJ3ZlLoJl+t6pk//\nnMceG5xlqQYNGrB48Vx69bqblJTKJCTcgNl8HpNpNh06xDBz5qw8fbCFhORlwvWsuDCbQ7Jc26BB\nA5KTD+DpyOOPLwCHcbvjqFmzph/2pfmTzWZj1qypHDp0iA8/nMxvv32LyWTiuuuaMHDgbsqUKRPs\nELUg0slYK1ATJ35BQsLAHMslJj7ExImvZZuMAdq2bcuxY/tZtGgRW7duxWYrSffua6hdu3aeYyxV\nqiRwKs/bZxQWdpqSJbNOjBaLhVtu6cG3336GyFP5ri8k5BP69Omrr64KsapVq/LmmyODHYZWyChP\n83YAK1BKAl2HVnTUrduS3bvfB3JqPt5NuXK3cPz4noII61/GjXufl15ai9P5VT73lIZh1ODnnxfR\nsGHDLEv98ssv3HjjXTidvwPh+ajvDDbb1axbt4wGDRrkYz+apuWHUgoRydU3Yn3PWCtQERER+HbV\neQqHIz+JKe/6978Pt3sJcDyfe/qO2rVrZJuIAVq1akX37jdisTyfj7oEm20I99/fVyfiK0BiYiIn\nT54kNTU12KFofqKTsVag7ruvB4aR/UAYAFbrNO6+u2cBRHS5qKgo7rjjDkJCJudrP+HhE3jhheyb\n2S+aOHEMpUotw2Qam6e6QkP/j/Llt/Luu2/kaXut8EtOTmbatGlcfXUroqJKU63aVdjtkXTpcjvL\nly9Ht0AWcbntfp3bBf1ok5bB+fPnxTBKCvyazeM5v4vNVkKOHj0atDi3bdsmhlFGYFueHi9SapqU\nK1dDkpOTfa7z4MGDUrFiLQkLe0J8GyJRBC6I1fqgVK9+VVDPV0Fzu92yatUque22flKuXC0pXbq6\nNG/eQb744otcnfOi4vjx41K3blOx22/0PnZ3cejQOIEJYrfXkT59HpC0tLRgh6pJ3h5t0slYK3Dz\n588Xm62seOayzTgesUtgodhsFeWzz6YHO0yZPv0LsdmqSO5nFfpWwsPLyvbt23Nd55kzZ6RXr3vE\nbo8W+DCLZ1JF4Lwo9b4YRhW5996HchwDuziJi4uT9u27it1eW5R6T2CnwF6BOeJwdJSyZavJ1q1b\ngx2m3yQkJEjt2o0lNHSYZD1WdLyEhraTsLCSEhpqlcjI8jJw4GM+DaCj+Z9OxsXEqVOn5K233pG6\ndVtIuXK1pHbtZjJixOty7NixYIfmN8uXL5f69VuKYVQTwxgghjFQ7PbaUqNGQ1mwYEGww/vbpEmT\nxWa7OMNOdqMnea5STaY3JCqqgqxbty5f9S5fvlw6deopYWHhEhHRRqzWRyQ09FmxWh+WiIhrJSws\nXLp1u0vWrFnjpyMtGlJTU6VVqw5isQzI8veh1BcSGVle9u3bF+xw/WL8+AliGLdkk4j/+fuDsgKb\nBPZKaOjLYrOVkf/9b2KwD+GKk5dkrHtTFzLTp3/BoEFPAN1JSrofqACcwmqdDszm3Xff5LHHHglq\njP60ceNGNm/ejNvt5pprrqFVq1YBeSzH6XQyc+ZMxo37lL/+OohSJqKja/H00wPp1atXtuMcL1++\nnKeeeoW9e/8kJeVhXK4BQEVA4RnMfxtW64eIfEWHDh0ZP34UNWrU8Evc58+fZ9OmTWzfvp2kpCTs\ndjsNGjSgSZMm3s5wV5aZM2fy4IPjSExcBWT9/LbJ9Bq9eu1j9uzPCi64ABARoqMbcvDgOKCDD1sM\nA5zAe96f92OzxfDll+/To0ePgMWp/VteelPrK+NCZPbsOWIYlQR+z+Kb714xjBryySdTgh1qoZeW\nliZz586Vzp17SdWq9cRkcojFcpPAd+KZrnC/wGwJD+8gUVEVZPny5Tnu87fffpP77hskNlukmExm\nCQuLFKVCpFSpqvLyyyPlyJEjBXBkV7bGjduKb/MAnxarNerv8Z1dLpcsXrxYOnbsIWXL1pQyZWpI\n69adZc6cOZKamhrko8raX3/9JVZrGR+uii8uWwUuna5xkdSq1UQPs1mACFQzNTAFOAFsy/BeSTzT\nkOwGlgJRWWxbIAdf1KWlpUmJEpUE1ubwz/a7OBylxel0BjvkQmvv3r1SpUpdCQ9vIzBYPBO5r8/m\nnP4ohlFGli1b5nMdycnJcvbsWd1hpgAlJyeLyWQWSPUpMUVE3CzffPONHD16VOrXby4OR2OBSQK7\nBfYIzJDw8HZSsWIt2blzZ7APL1M7d+6U8PDaPiZiEfhLoMIl77nEbq8p69evD/bhXDHykox9fbTp\nU6DLJe8NBX4QkTp45gAbmqtLcu1fvvvuO9LTqwGtcyh5NdCSr77K74AUxdOZM2do06YjR44MIT7+\na+ArYAnQPJutOuB0zqZnz74kJib6VI/FYqFEiRKEhupB7ApKUlISoaFWwOxTeZFwTp06RatWHdiz\npzsJCZuAgUBtPHMI9yE+fiXHjr1ImzYdOXToEBs2bKBPnwE4HKUJDQ0jIqIcAwYMZtu2bYE7sGxE\nRUWRmnoa34dnPQ6UuOQ9EyZTC3bv3u3f4DS/8ikZi8gqPBNuZtQduHhD5jNA35DIh++/X058fG+f\nyiYk9Gb+/OUBjqhoGj/+Qy5ciMHtHgx8AtwONPJhy/a43a354osvAxuglmfh4eF4ktJJn8ortZ/l\ny3/ixIlWpKe/jOce/+VEBnDhwgBiYrrRvn0vZs+uS2Lib7hc8cTH/8r06eW49trODB36ysXWvgJT\nvnx5atWqCyz0cYvpQK9M3k/VXxwLufwM+lFORE54X5/g39PvaLkUH5+E70MhOkhMTApkOEWSiPD+\n+x+RkjIEEOBj4FGft09MHMx///tRoMLT8ikkJIReve7CZJriQ+kthIQcZdGiZaSk5Dyymdv9JAcO\n7MPpXInb/QJQGbAA1UlPf5WkpC188MF83nhjVD6PIvdeeGEwdvtocr46Pg58Djx8yfuJpKXF5mkG\nM63g+OWrkoiIUirLr4wjRoz4+3VMTAwxMTH+qLZYqVGjIqGhu0hPz7msybSbGjUqBj6oIiYpKYnz\n508CTYAE4DTQOBd7aM3Bg38EJDbNP5577nHmz++G03kPnoSZGRc224v07HkLc+ZsAOr7sOdSQEdg\nHZBZT/iyOJ0LeOONa3jssYcpUeLSpuDA6dOnDxMmTGXz5odISfmYzD+2TwI3A0OAKv9ao9Qkrruu\nLVWrVg18sFeo2NhYYmNj87cTX28uA9X5dweuP4Dy3tcVgD+y2C6Q98mLjV27dnl7TSbn0EEjXQyj\nmu6MkYmUlBQxmULF8/zpOYHwXHR8EQGnhIZagn0YWg7eemu0GEatLDrlHRWbrZdcf31nmT9/vkRG\ndsrF7/9BgY+yLWMYd8vo0f8t8GOOi4uTtm27iMNxlSj1P4HjAk5vZ7RnBUoJDJfLe13PlPDwsnrw\njwJGADtwZWY+0N/7uj/wbT72dcWrU6cOLVu2wGx+PdtyISFjqVWrEs2bZ9ch6coUFhbG1Ve3BBYA\nF5/BPZaLPeymZEnd4lDYXLhwgVmzZvHMMy/QqVNvFi9exVVXVcYwOmEY1wAjgTex2+/Ear2K++6r\nwrJl8ylbtixud24m+zjG5Z2f/s3p7MV3363Ix9HkTXh4OCtXLuK778bTpUsshlGfkJAooqJi6Nfv\nDJUrl8ThmAe8D8wFJhAe3pIKFV7ip5+WUL++L60DWjD5NOiHUmoG0B4ojef+8CvAPGAWUBU4CNwp\nIucz2VZ8qUODkydP0qxZW06c6Epa2jAg42Tj5wgJeZeSJb9gw4afdJNTFr744gsGDfqAxMSf8DTZ\nVQJe9mlbi+UxnnmmNG+8oeeaLQyOHj3K0KEjmT17Fmbz9SQkXItIXTy9qU9htW5E5HvsdhPXXdeI\nm2/uQt++fYmMjAQgNTWVqKhKJCUtAZrmUNsx4Co8H2WR2ZRbSvPmo1i/fln+D9CPRITly5czZcoM\njh07TZkyJbjvvtvp0qULISFZD46iBYYe9KMYOHXqlNx11wNitUaJw9FLLJbHxeG4U6zWKOnRo58e\nWCIH6enp0rlzD7HZugksE6gocNaHJsrDYrWWkL/++ivYh6CJyPTpn4vDUUZCQ18UOJbN780znrnd\nfrV07tzj70E+REQGDBgsoaE1BHpm0nx76TJEYJAPfyfjpWfPe4J4ZrSigDw0U+vhMAups2fPsnDh\nQs6ePUtkZCQ333wzZcuWDXZYRUJqaipPPfUiU6ZMIT29POnp4XjHpclii+MYRmdefPFuhg9/oQAj\n1TLzxhujePPNj3A6v8b3DngpmM0vUqHC9/z66wpOnTrFtdfeRFLSRuA2oBWeISIvvUoUTKa3ERmF\nyC4gu/8xweFoytdfv0Pnzp1zfVyX7U2EtWvXsm/fPkSE6OhorrvuOkwmPbNtUZeXK2OdjLVi68KF\nC3zzzTd8+OGnbN68n7S0pxB5gH/uC54gJOQTLJbxPP30I7z22ssBGRf7SpSamsqmTZtIS0ujYcOG\nfzcd52TGjK8YOHAYTucqPON/505o6Ehq1fqO669vxaeflsblGgGcB/oCO4BBeO64mYBfsdkmUL16\nJEeOHCIubjqQXZL9ksqVR/LnnzvzlTBdLhcTJkxk1KjxnD9v4p8m9M1ERKTx7LODeeKJx/RzwUWY\nbqbWtCysWbNGeva8WyyWcHE4aojDUV2s1ki5++6BsmnTpgKNxe12y7lz5+Tw4cOSkJBQoHUHmtvt\nlrff/q9ERpaXiIjGEhnZWqzWKHnooSckKSkp222PHj0q4eFlBTbkshd8xsUthtFVKleuL555fzOu\n2yjwkEBrgZYSElJPhg4d+vfcyIZRWmCGt+k743ZpotTHEh5eVrZs2ZKv85OSkiIdO3YXw2gnsPKS\n5nO3wCoxjBukffubczxfWuFFHpqpdTLWrihxcXGyZ88e2bdvnyQmJhZo3efPn5dx496XypXridns\nEJutvISGWqVp0/Yyc+bMQj1hga+efXaYGEZTge0ZksxRsdnukPbtu4rL5cpy24ceelzM5mfzkYgv\nLn+JyWQXGJNtufDwLjJnzpy/6//111+lTp2mYhg1RKlhAmMlJOR5MYxK0qjRdfL777/n+/z07TtA\nbLbbBFKyiS1VbLbbpXfve/NdnxYceUnGupla0wrAmjVruPnm3qSnt8fpHAy0wzM8YyrwLQ7HeMqW\nPUds7EKqVKmS/c4KqaNHjxIdfTUpKbv595MAAOk4HM2YM2cUN91002XbJiQkULZsVZKStpL1YB6+\nCwu7F5HVpKUdyKLEYWy2hpw48ad3mE0PEWHdunV8++13nDlzgXLlSnL77b1o1MiXIVWzt2/fPho0\naEVy8kHAnkPpJGy2GmzaFEu9evXyXbdWsHQztaYVQuvXr/c2gS7Otnk1JGSUVKhQU06ePBnskPNk\n7NixYrUOyFNP5Pnz50tEREc/XBVfXJaJ2VxWQkLezmRdvBjGDfLMMy8V6Pl54olnxGx+zudjCA19\nSQYNGlKgMWr+QR6ujHW3vSuciHDs2DEOHjxISkpKsMMpdkSEXr3uxemcCFx+RfgPhcv1HKdP9+Dx\nx58rqPD86vTpMyQnZ3dVX5WTJ89muubXXzeQmNjSj9E0BxKpXHka4eFt8Uw8t4iQkDcwjKvo2TOa\nd955zY/15WzevCWkpd3lc/n09Lv47rslAYxIK0x0Mr5CiQhTpnxKzZqNiY6+hgYN2lOqVCUef/xp\nTp06Fezwio0VK1Zw7pwZ8G1GrrS0ocyfP48zZ84ENrAAqFu3Dg7HuizXh4b+SoMGtTNdt337flyu\nOn6MJpLQ0HBiYxcydepT3HzzYlq3/oD+/Y+yevU8Pv/8kwIfDCMxMQ7PGNi+KuXdRrsS6GR8BRIR\nBgwYzBNP/I8DB94lOfkkiYl/kpi4no8/TqVRo9YcPXo02GEWC++9N4mEhEfJavq+y5VGqVuZPv3z\nQIYVEL179wbWA79ksvYYZvMkhgy5dEYhj/R0F36at+ZvJlMoJpOJXr16sXDhTNau/Z7Jk8fTpEkT\nvx7c2iIAACAASURBVNbjK7s9AsjNl6yzGIavM7lpRZ1OxleguXPnMnv2apzOWKAT//wZ1CAt7X+c\nOnUv99zzSPACLEZ27twNtMjVNklJLfj99z2BCSiAbDYbM2Z8imF0R6n38MwkFA98gd3elhdf/A+V\nKlXigw/+x113DeDWW/sxcODjLF26lPLlS+IZaddf0klNPVugsyvlpHv3mzCbZ/lcPjR0Jt27dwlg\nRFphop8qvwK99db/SEx8mazmT05Pf46ff67KgQMHqFEjs+nkNF+5XC4uH/UpJyGkpfkwl2Yh1K1b\nN1atWszrr/+XxYtH4HKl0bTpdTz33DssXLic8uWrYzJ1xunsCBjAMWbOfI6QkONYLDH4r9vCH5Qu\nXelfPaWD7cknH+GTT9qQlvYqnmPPTjJm8yf85z+xma49efIkU6dOY9s2z5e2xo3r0b//vZQuXdqv\nMWsFRyfjK4zb7Wbz5lVAdh1DDEJDO7Nq1SqdjPOpUqWKHDy4G88cy74xm/dQvXrRnT2qadOmfPPN\nF3//fOHCBdq06cT+/fVITt6BZ8bVfyQkPAPMxDM6VgpgyXcMJtNCbrihXb7340+1atWie/db+O67\ne0hKmolnwovMpGOz9adLl46XPdaUkJDAwIFDmDfvG6AXycktAOHrr39m2LCR3HVXHz788D1sNluA\nj0bzN91M/f/t3Xd4FMUbB/Dv3OXa3l0KCQkJNXQUkSaCVOlFaSqCoAiIgPJTREWKgmKhCSiiNEFE\nLIgoRpqAGLqANIN06U0wBEi7lLvv7487IEDK1VwS5vM8+0Dudnfe2WzuvdmdnbnLXO9Gn1drjQyA\nzWbLn6CKsBdffBpm8+cubJEKtXohnnmmp89iym+dO/fEP//UhcXyJW5PxHYCQHfYx6Fe4oUSrTAY\nZmLIkAFe2Jd3LVgwEw0aZEBRWgPYBCDrGAwEsAWK0hZ1617DN9/cet4kJyejQYOW+PlnwGI5Botl\nLoCBAAYhNXUB0tKOYvHia2jcuC0sFku+1UnyDpmMC6iUlBTs378fe/fuxYULrszJmju1Wo2KFWsC\nyG0KuAyQ61CnTh2vlXu36tq1K4SIA7DbyS2+RJ06dVGhQgVfhpVvdu7cie3b45CWNg15d2J7C8Bw\n2O8zu0+tnoLq1cvjgQdcu1efH3Q6HX799Sd88EFXREX1g8lUA0ZjbxiNvWE210RkZG+8++4j+O23\nGOj1+lu2ffHF13D0aBVHEs7uXngYUlMXYv/+CAwb5ty0oVLBIUfgKmAOHDiAyZM/xTfffIOAgHAI\noUVa2hnUrFkHb7zxAjp27OjxIxlz5nyOIUO+QkrKb8juToUQc1CjxgLs2bPRo3Iku6+//hbPPz8C\nKSnrAJTPZc21MBqfwubNa7wy4lNB0LPnc1i0qAKs1hFObtEXQAqAb+F8D/SsdkFRWmPfvh0F/hYL\nSWzatOmWWZsaN26c7SQUCQkJiIoqD4vlIICIPPZ8EopSG+fOHcMff/yBM2fOQK/Xo3HjxnIe9Hwi\nR+Aq5D7/fB4NhuIMCBhN4HSW0XgsBL6m0ViHLVo86vGYyunp6WzatD0Nho4EjmYpJ4lCTKXZHO6V\ncXilmz755DMqSgkKMZlA/G2jLR2jRvM6TabiXL9+vb9D9arAwAgCJ1wYOesagXJUqQYSyHRx1K0/\naTCU4A8/LPF3tb1uxowZVJQnnT4WWm1rKkoIzeb6NBr7OOZEL8bmzTvKv+18ADdG4JLJuIBYtOh7\nKkopAgdz+SNLo17fnS1bdsx1wH1nWCwWvvrqCJpMYQwMrMegoObU60P58MOPcv/+/V6qlZTV9u3b\nHTNHBTEoqAUDAzszMPAhGo2hHDz4VR4/fpz//vsv3333A1at+iDDwyswLKwcy5evxUGDhvDQoUP+\nroLL1GotgRQXk+pUhoVF02hsQOCAE+tnUK0eR0UJ5ZIlP/q7yj7xyiuvExjnwjEcSWDQba8lEZhG\nk6k4//jjD39XqUhzJxnLy9QFQHp6OsLDy+Lq1Rjk/UxqOozGeli06H106NDB47JTU1OxY8cOpKWl\noWrVqoV2koLC5OLFi9i9ezeSk5MRHByM+vXr4/z583jttdFYtWoFgMdgsfQCUBL2jnb/QaP5CQEB\n83D//fdj0qTRaNSokX8r4SRFCUFq6hEArjxyMwUDB55E1aqVMHLkGAjRBMnJzwB4EDfnOLYAiINa\nvRI63RzUrn0vFi6chbJly3q7CgXCsGEjMGmSCcAoJ7cYAftkFG9m894yBAc/j7Nnj0JR8nrESnKH\nvExdSH333Xc0mZq58K33CzZu3N7fYUtesmXLFgYGRlClep/A5Vx+7xYCX9JgCOf8+Qu8Vv6VK1c4\nduz7jIqqTL0+kNHRNfjxx9O8Mp/uQw+1IfCVSy1js7klv/nmG5JkYmIiZ86cxfr1W9NoDKVWa6Ze\nX4wBAXqWK1eD/fq9yD179ngcZ0G3cOFCmkxtXDiOjQgszfF9k6kD582b5+9qFVmQLePCqUWLzli3\n7nEAvZzcIhVabQQuXDjp9RGGrl27hkWLFuHQoaPQ6bR4+OGmaNGiBYRwpzNNwWGz2XD48GFcvXoV\nRqMRlStXhlar9XdYiIuLw0MPtUBS0nwA7Z3caj8MhlZYsGAaHn/cuTGvcxIfH48HHmiG8+erw2IZ\nCqAigL+gKONw771pWL9+hUfPrMbExKBHj/eRkrLNyS0OITCwKS5ePAmd7tbnjUni2rVryMzMhMlk\nuuP9osxisaB48TJIStoC++8oN38B6ADgOHIeSiIGNWtOxu7d670ZpuQgW8aFVLVq9Qlscqn1YDJV\n9Oo9RKvVymHD3qLBEEKj8XECHxAYTZPpPpYsWZnr1q3zWln56fLly5wwYRJLlKhAo7Ecg4Lq0myu\nwsDACA4bNoqnT5/2W2yZmZmMiqpI4GsX76mSwE4aDMV49uzZO/YbHx/PSZMms2rVegwLK8eoqKp8\n5JEn+fvvv9Nms92ybs+ez1GjGUzAdtv+rdTrO3P06LEe17FYsVJOto4zaDC058iRYzwqMy0tjcuW\nLePs2bP51Vdf8dSpUx7tr6B4882xVJSmBFJzOYZJBOoTmJzHsT7E8PAK/q5SkQU3WsYyGRcA99/f\nlMBvLn0YG42leeLECa+Ub7PZ+PTT/akojQmcua0sG4EVVJTiXLNmjVfKyy9///03ixcvS4PhKQJb\nb0s4+6nV/o8mU3GuXbvWL/H98ssvNJvruZGI7YtON5Bvvvn2jf3ZbDa+++546vXBVJRejnPqKIF9\nFOITmkzVWLHi/Txy5AhJ8tq1a9TrgwhcyKGMvQwJKXlHAnfVkiVLCBgJzM8m6V9fkqnVdmLjxm2Z\nnp7uVjkZGRkcPfpdBgWVoNnciIrSjyZTN+p0IWzZsjMPHjzoUT38LTMzk5069aCiNCLwZzZ/p1sJ\n1CDQL5fjfH3ZxTJlqvu7SkWWTMaF1MCBLzMgYJQLH8QHaTaHMy0tzSvlr127lkZjZQKJuZS5lqGh\npZmRkeGVMn3t1KlTLFasJIWYn8ex/J2KEsZt27ble4yNGrV1JCj3kjGwl8HBUTeS16uvjqDRWJPA\nqRzWt1Gl+ozBwZE8evQo9+3bR7O5aq5lBAQovHbtmsd1LVfuXgJlCNQkMJvAcceXgD0EXiFg4mOP\n9aLFYnFr/xkZGWzdujMVpRWBv2+rRyKFmEyzObzQ31+2Wq2cNGkKw8LK0mx+gHr9C9TrB9FsrsXi\nxaMZEGAkcDXPcycg4E327z/Y39UpsmQyLqT2799PgyGCQJpTH8Ja7RC+9toIr5Xfpk1XAjPyLNds\nbsClS5d6rVySPH/+PKdMmcKXXnqVo0eP4fbt272y3+7d+1CtHulkUvuK1as38Eq5zjp9+jT1+lC6\n/tjP7b+ThoyJieGqVauoKOUJXMpzG5VqOqtUqcMzZ85QpyuWy3n3HzUaxStfwNavX0+DIYzAGAId\nCZQmEEYgmgEBoRw79r1ct09ISODq1av566+/Mj4+/o73x479gIrSkkB6LnX/nuHh5dxueRckmZmZ\nXLlyJadNm8ZPPvmEa9asodVqZZs2XR3Psud2DlyhwVBCPsLoQzIZF2KNGrWlRuNM8thJRQnz2iVq\nktTpTAT+c6Lsqezb9wWvlJmcnMzu3ftSpwumXt+PwASqVCNoNEazatW6jIuLc3vf8fHxjsuv/zqZ\n1DKoKKW4a9cur9TNGRs2bGBQUEOPErH9i9n/+NFHH7Fx4/YEvnByOyuNxkrcvHkza9duypzu56pU\nH/Cxx572Wp1///13li9fg0ZjRQYGdqbZXIthYWVy7RmenJzMvn1fpF4fzMDApgwKeph6fTCfeuq5\nGy329PR0hoSUJLA3x78Z+6XbMAKCGo2RDRq05tKlS5mZmem1+hUEBw4coNkcTuDHHBOxojTh88//\nz9+hFmkyGRdiFy9eZKlSlanRDGf2HTRsBNbQYIjgkiXeHWFIpQrIpXWUdZnDJ5/s63F5FouF9eo9\nTL2+F4Ert5WRSSHmeTQK2Jw5c6go3VxKairVaA4ePNTjujlr+fLlDApq63EyFmIkhw4dSr0+jK60\nsoWYwi5denHTpk1UlOIEfiJgvfHlBJhLsznc6/dZbTYbd+zYwSVLljA2NjbXZJiens769VtQr+9B\n4HyW+C9Sp+vLGjUaMDU1latXr6bZ/GA29Uwm8Djtl8c/oP3yfTrtj48toNncgGXKVL1xD72o+PPP\nPxkaWsrxKNSPBPYT2MGAgJE0GCL4/PP/K3JfQgoad5KxnEKxgChevDh27dqEJ5/si61bS8Nq7YOM\njKYAtAD+gck0B4qSgAULvkSbNm28XHZZ/PtvHIDcJ4bQauNQubLnY9vOnDkL+/ZpYbHMx52zR6lB\n9kFSUgZ69hzo1vjYZ8+eRWpqZZe2sdkq4/jxZS6X5S6TyQRPJ0QAAI0mEYmJFuh0dWGxOP8IEtkU\nu3d/iYYNG2L58u/x3HNDcOHCqwgIqID09L9RqVI0FixYjSpVqngcY1ZCCNStWxd169bNc90lS5Yg\nLi4VFstXuPU8KY60tDk4erQtFi5cCK1WC/L2x30yAHQGUBzAYdw6LWMIgKeRmPg0kpNn4sEHm2H3\n7i1FZtzmOnXq4MyZI1i8eDGmTZuOc+fOQqvVoU2bZhgyJPaOaRmlAsLV7O3qAtkydtmRI0f48suv\ns0GDtqxTpzk7d+7JFStWeDwEZk7Gjn3fcak4t9ZUMvX6UI8vj9tsNkZFVSYwgcALBJ4i8DyB725r\nnafTYIji3r17XS7j/fffp0o13MVW5nx27tzLo7q54tSpU166Z9yIw4cPZ2BgRxe33ceSJavdiMdm\ns3HPnj1ctWoVDxw4kG/HITf16rUk8H0udVjOe+6pz8WLF9NsfvS29z4k0NrRys/9WKjV49mwYRt/\nV1cqQuCPlrEQ4gSAawCsADJI1vN0n3e7ihUr4qOPJuZbeQMGPIfJk++HxRIDoGM2a2RCr38e7du3\n83i4walTP8a5c+cB/ACgB+zDJF4FMAvAEADDHP9qYLU+htWrV6NGjRq57vPkyZOIiYlBQkICgoOD\nERoaCqPxVyS60PDU67eiVq38azGULl0atWvXxZYtiwE84+Ze4qBWH0O7du9j+vTlLm57HOHh4Td+\nEkLg/vvvL1CzRZ06dRLAfbmsUQPnzp1C06ZNkZHRH0AC7K1eG4DPAHyDnAe9uMlqfRk7d07GkSNH\nUKlSJS9ELklucDV7377APsxLsVze9+1XEMkrtm3bxsDACGq1L/HmZBWZBGJoNDZmgwYtPZ4t6r33\nJlCvj6b9ecjsW2tALQL2QSiEGMmxY3MedOLUqVNs3rwj9fpijpb9m9Trn6NOF0K1OtCxP2daiVeo\n0wVnO4CGL8XExORwr9O5RacbxFGjxnDPnj0Uwkxgh9PbGo2dOWPGzHytr6vuvbcBgVW51GMTy5W7\njyTZtWsvqtVvO15f6ziP8nrW9uai0QzjkCHD/FxjqaiAGy1jbyXj0Fze93G1JW85c+YMX399JAMD\nw6nRGKlWa1mtWj3Onz/f48dBVqxYQUUpQ+BsHh+MVwncT2AWTaZO/OKLL7Ld3/HjxxkaWppq9bu0\njzqUdR9JBFoQaMWbnZJy+yB+nY888qRH9XNHZmYmIyMrEPjWjWS8mwZDMZ4+fZqVKtWivaNSdye3\nPUyDIYSJiYn5XmdXTJ/+KRWlc471MBh6cdy4CSTJEydOMCQkikLMIfAZgf4uHs/v2bJlVz/XWCoq\n/JWMjwHYDeBPAP2zed/nFZe8y2az8erVq16ZKOC6unUfdiHpbCJQnnp9cI4DTtSu3YRq9Ye57MNC\noAqFeJI59xS3MSBgLKOiKvLixYteq6sr9u7dS6OxOIGVLiSOAzQYSnLRou+5fft2Go0VaO8hfA/z\nnmbvLBWlCj/55DO/1NcViYmJLFOmKgMC3uGt934zqVJ9yIiIaF6+fPnG+ocOHWLZstWo05Ui8JyL\nyXgJH364sx9rKxUl7iRjb/SmbkjyvBCiOIA1QoiDJG/pAvv222/f+H+zZs3QrFkzLxQr+YoQAoGB\ngV7b34EDB/D33/sBdHVyi4cAqNG+fRuYzeY73t27dy8OHjwGq/XlXPahA7ANQpSHTheNjIyBsFp7\nwN679hqE+BlG42eIjg7Cr79uQPHixV2tllfUqFEDq1cvRbt2XZGY+ArIAQCCc1g7HcBiGAxD8emn\nE9Gt2xP47rvvoFLVgf1e6SoA7QDsAvAKgPoAro9VfxVCfAWDYQJGjBiMwYMH+bhmnjOZTNi69Td0\n7PgUDhyYA6u1E0gVNJpfUL58JGJi1t0yUUrlypVx7Ng+vPvuu3jvveXIzHS+LJXqACpUKOWDWkh3\ng9jYWMTGxnq2E1ezd24LgDEAXr3tNZ9+A5EKvvnz59No7OliS+UNvvPOO9nub+jQN6hWv+nUftTq\nN9mr17Ps3XsAQ0JKUaczMzCwBDt06Mb169d7PO6ytxw+fJiPPtrdMQhKfwIbaB9X+jizPiNar14L\nrl+//sZ269ato9lcO0udr9I+SUBFAlUJNKZe34A6XTDbt3+CGzZs8GMt3bd7925OmTKFkydPznOU\ntrS0NAYGRtD+fK0z51omFaUs//zzz3yqjVTUIb9bxkIIBYCaZKIQwgigNYB3PPt6IBU1qampsNmM\nLm5lwhdffI2VK7cgKioczz3XA23atIFKpcLp0//Cam3k1F6s1mikpZ3B99/PdD3wfFSpUiXExHyL\nCxcuYPbsuVi4cCiuXImH1ZqJwMAQtG7dFK+8cuczok2aNIHBcBmJib8BaAEgEMBQ2Huk/wKN5il8\n++3XqF+/PkqUKJH/FfOSmjVrombNmk6tq9VqMXDgc/j44/FIS5uPm1cHcvItypYNR506uT9nL0m+\n5Oll6ggAPznmug0A8DXJ1R5HJRUpoaGhCAg46+JWJ3DiRG2cONETwEmsXv0mgoJexrJl3yM42ATg\nipP7SUBQkKtfBPynRIkSGD16FEaPHuXU+mq1Gl99NQtduvRASso7sM+JrQWwFIryBiZOnIxHHnkE\ny5Ytw4IFS3DpUgJCQ4PQq1cXdOrUCRqNxpfV8ZuRI1/H4sWNcPr028jMfBs5J+Q1UJRX8NVXq/Ix\nOknKhqtNaVcXyMvUd70rV65Qrw8mcM7Jy4YptI8j/M8tna2A72gyFeenn35Kk6kO8350xUaTqQ6X\nL1/u70Pgc9u2bWOrVl2oVmuoUqlZp04zLl++nOvXr2dYWGmazQ8RmEUghsAcms2NWaxYyUI3LaYr\nzp8/z6pV69BkakZgSZZOYDYC22gw9KbZHF5oL91LBRfcuEwt7Nv5jhCCvi5DKvieeWYAvvmmOKzW\n95xYezaApQBW3PGOELNRs+ZCnD17BhcvzgLQKpf9rEV4+PM4d+4I1Orbh90smq7/YatUKmzatAlt\n2nRBSspCANkNofo79PrH0LdvD5QqVQplypRB586dYTQWnisJeUlPT8ePP/6ICRM+w/79e6DVhiEz\nMwmBgSYMHToI/fr1QVhYmL/DlIoYIQRI5nV/5FauZm9XF8iWsUT7c8FBQSUI/JBHa3aDo1W8M4f3\n02kwRDomgwgnsDGH9TbRYAjn6tWr/V11v7BarYyKquRoCWd3fPYTaEYggkA/qtXDaDJ1oKIU47Bh\nbxW6iQRsNhtjY2PZoUM3mkzFqVIFUKNRGBFRgcOHv8XTp08zPj6eR48e5blz53w2tKwkke61jGUy\nlvLNzp07GRwcRa32f7w5ytf15RSBEQSKE1ida8IOCHiZ48aN48qVKxkYGE6TqTWBrwn8RuBrmkxt\nGBgYzpUrV/q7yn6zevVqmkw1criUv8+RhD/jnfP/nqCiNGfXrj0LTcJas2YNy5S5hyZTNQrxCYEz\ntD9bnkhgN3W6F6nThbBdu8d56dIlf4cr3QXcScbyMrWUr86dO4epU6dj1qy5EKI07I+nX8W1a3tB\ntgMwHsDtM/Dc7j0MH56KcePeh8ViwQ8//ID585cgPj4BxYoF49lnH8MTTzwBvV7v+woVUAMGvITZ\ns0sDeD2bd+sBGACgXw5bp8JobIh584ajW7duPovRGxYsWIiBA19DaupcAO2Rc0etRGg0YxEe/jO2\nbFnr8xmaSGLjxo2YP/87nDlzEXq9DvXqVUf//n0RERHh07Il/5OXqaVCw2KxcOPGjZwxYwY/+eQT\n1q7diDlNcn/7otP155QpU/ItzunTP+U999RnZGRldujQjZs3b86Xsj3RteszBOZlc/y2E4imfdzx\n3I7zYtau3dTf1cjVqlWraDC48jwxqVZPYrly9zIhIcGncZUtey9NpqoUYiKBRQS+ol7/HPX6YHbt\n2sun5Uv+B3mZWioMMjMzOWPGTJYrV52KUoZBQQ2o05WjECEEJtE+lGVOH6hJ1OuL8eTJkz6P02Kx\n8MEHm1NRWjsunR+gENOpKJGcN2++z8v3xIsvvkKVamw2x2+043ZAXokrnRqNwqtXr/q7KtnKzMxk\neHg5x60J5xLxzS9zz3D48Ld8EtfChd84viAsz+EWQQK12sGMjr6X//33n1fK3LdvH6dOncqJEydy\n7dq1heb2QlEmk7FU4GVkZLBdu8eoKA1pn13n+geWjcAfBNoQaMo7J3+g437xaLZqlT9jCE+b9gkV\npQ3vnGziAA2G4ALdutm6dSuNxvLZxD6UwESnkpbBEMFz5875uyrZWrZsGc3mB1xOxPblbwYFRXo8\n+cntdu/e7ehUGHfjiyNwhMBh2u9f34xBo3mdDRq09Ki8c+fO8cEHW9BgiKRON4gazSs0me5jyZKV\nC8XVm6JMJmOpwBsy5A1Hgsup9Wsl0NOxZH09lQEB7zAiIjrfEkTVqvUIrMk2TqPxcc6dOzdf4nCH\nzWZj5cq1KcTnt8X+IZ2b0eg/arUmr04W4k2NG7cj8IWbyZg0m5tx0aJFXo3piSd6U6WaSGAXgb4E\nggiUJ1CBQCCBpx1fOEn7EJxl3B6CMyEhgaVLV2FAwNu8tROejcBSKkoYd+7c6dX6Sc5zJxmrvHjP\nWpJylZiYiFmzZiMl5XPYJ3LIjgrATNifMx4CYCo0mpeh15fGQw9tw65dmxAZGZkv8V6+HA8g+44+\n6ellEB8fny9xuEMIgZ9++gpm8ygIMQvA9VkTegD4AcC1XLdXqb5Ax45dC2QnOJvNhs2b1wBwv3NZ\nYuKT+PFH7426FR8fj5iYpbDZrgLoAHsnxMMA/gFwFPbJ7e4H8ASAkQBUSEsbiMmTP3OrvM8+m4VL\nl+ogM3MMgKyjqAkAnZCS8i6GDHnL/QpJ+U4m4wImMzMTS5cuRcOGbREUFAmNxgBFCUZUVBUMH/4W\nTp8+7e8Q3bZ48WKoVE0B5DU7jgkaTW889NAeDBhwAqNHhyMubivWr1+OqKiofIjUrlatmgDWZfMO\nodP95vRYyXlJT0+HzWbzyr6yuueee/DHH7/j3nvnQ1HKQ6UaCWAJgGjYe1PnVOZB6PUfYsSI3GbF\n8p+kpCSo1XoAigd7CcfFi5e9FRK2b98OIBLAT7DPJjsCQHiWNUIBvApgJ4BfAbwLq7UL1q1b71Z5\n06fPhcWS2++nN3bs2IazZ10dhlbyF5mMCwir1Ypx4yYhIiIazzwzCVu2PI1r13YiMzMeqanHcP78\nd/joo6uoXLkmWrXqgv379/s7ZJcdPXoMycm1nFo3I6MOoqLKYebMj/Hmm6NQsWJejzt536hRL0NR\n3gNwJMurhEr1EYoXt6JFixZu7zs5ORnvvz8e4eHRMBiM0GoNeOSRJ7Fjxw6P486qWrVqiIvbio0b\nl2LYMBX69DmIIUOao2rVwzAYOgLYCoCOtRMBzIDB0AyffjoRtWvX9mos3qJWq0FaPdyLFQEB3phB\n1u7s2bNISzsN+6hxuX1hLO5Y52MAiUhNTXKrvIsXTwCokcsaBmi1ZTFt2jRUrlwHer0ZWq2C0qWr\nYfz4ifjvv//cKlfyIVeva7u6QN4zzlNKSgrbtOlCRWlKYHce97sSKcRHNJmKc926df4O3SVjxrxN\nYJST9/VmsGfP/v4OmbNmfU69PphG45MMCHiNZnMtRkdX5/Hjx93e59WrV3nPPQ9Qr3+MwA7Hfb4E\nCvExFSWcixf/4L0K5CA1NZWTJk1mZGRFKkopms3VqNMFsU2bx7hlyxafl+8Jm81GrVYhcNHte8ZC\nTGKfPoO8FlOfPv0IdHchhv8ReJ6RkZXdKs9sLk7gWC77X03AQIOho+P/V2nvRLaVBsOz1OuD+cUX\nX3qt/tKtIDtwFT6ZmZls06YL9frHmfsjPbcv66goYXnO7VqQrFq1iiZTTeY9wQNpMnUoMB2k4uPj\nOXv2bI4fP56//vqrx4+O9Ov3InW6vjkch900GELybaQoq9XKo0ePMi4ujvHx8flSpjd069abKtUE\nN5OxjSbTPfz999+9Fk+5cjVoH8rV2Rj2EIjgwIEvu1Ve//6DGRCQ05zesQRCHf/mVP7fNBhKqbd3\n9gAAEutJREFUy4TsIzIZF0Ljx39IRWniYiK+vvzMkJAopqSk+LsaTrFarSxRogKB9XnU6xCNxlAm\nJyf7O2Svu3btGg2GEAJnc6y/wfAMx4+f6O9QC7Rt27bRaHRm8JLsk1Xp0tVos9m8Fk9QUCSB0y7E\ncJWAlgcPHnSrvMOHD9NoDCPw+237zSBQksCvTsTwN/X6YK897yzd5E4ylveM/chqteLDD6chJWUq\ncu5dnJuOyMiogcWLF3s7NJ9QqVSYPn0iFKUngIM5rHUWitIJH3zwDhTFkw46BdPff/8NrbYCcruv\nmJraAWvX/pF/QRVCDzzwAMqUCYMQc13cMgOKMgbDhr0IxzzsXqHR6ACkubBFGtRqLapUqeJWeZUq\nVUJMzCKYTE9AUbrD3jFvBQICugIoBqC1E3u5B0J0xNy5X7gVg+RdMhn70fLly5GWFgXA/Y4ySUkv\nYPz4T70XlI899lhXTJ/+HgyGhtDpBgPYBeAigL+h0QyHXl8TI0Y8i5deetHPkfqGvfNRRh5rZSAg\n4O6Y8tFdQgj8+OMCmExvAfjZya0yoNf3QYMGgRg4cIBX46lVqxaA31zYYi1q1XrAozKbN2+OM2eO\nYty4hmjSZD7q15+GyMhjAF5zeh+pqQMwfbpMxgWCq01pVxfIy9Q5atSoHYEFbt73ur7YBw/YtWuX\nT2NNSkrivHnz+PrrwzlhwgSeOHHCo/2dOnWKw4e/xaioyjQaQxkeXp6DBw/l4cOHvRRxwWSxWGg2\nhxM4lOPv1GjsxBkzZvo71EJh+/btDAqKpEr1NoFLufyd/ElFaclmzdozKSnJ63HkPkvWnYvZ3IhL\nlizxehzly9ci8KcLnx9XqNOZvB7H3Q7ynnHhEhJSksBJD5MxaTL15Jdf+q4jxqJF31NRitFo7Ejg\nPWq1L1CvD2WfPi8wIyPDZ+UWVcOHj6bB0J72af5u/32uYFBQBBMTE/0dZqFx7Ngxdutm7yFsMDxN\n4Cfa57n+jcBsms31GBZWluPGTfTZ+Wq1WlmhQg2qVJ/k+fcqxJeMjKzgk1gqVqxD+2Qgzn5+XKZe\nH+T1OO527iRj7z1oJ7ksJeUqgGCP95OREYwrV654HlA2NmzYgD59XkJKym8A7INcpKcDwDgsWtQF\nJtNITJs20SdlF1Vjx76JXbu6YfPmhkhOfhVAAwDx0Gq/hE73HVas+Akmk8nfYTrlr7/+QmxsLFJT\nU1GiRAl07twZQUFB+RpDdHQ0Fi36AvHxH2Lu3C8QEzMHCQkJ0Gp1KFWqBAYMeAvt2rWDWu27S/8q\nlQpr1izFAw80wZUrybBah+DOfiAZUKlmwmz+AL/9ts6rzzlfV716Vfzzzx8gnb0Evg3lyrl331ry\nMlezt6sLZMs4RyZTGIF/PW4ZK0o/zp492ycxNmzYlsD8HMr+l3p9MC9fvuyTsouyzMxMLl68mA8+\n2IohIaVYsmRVDhs2iqdPn/Z3aE7ZtGkTa9RoSEUpSZ1uIAMCXqfR2Jl6fTD79BnEa9eu+TtEvzh5\n8iQbNGhFvT6cAQHDCSwm8APV6rdoMESxdu0mPHLkiM/Kj42NpclUzenL5SbTo/z88899Fs/dCm60\njIV9O98RQtDXZRRWZctWx6lT82Cf7N19gYGtsGDBYHTq1Mk7gTmkpKQgKCgUmZmXARiyXcds7oTZ\ns3uge/fuXi27sCKJTZs2Yc2atUhMTEWZMpHo3r17vo2nnR9WrlyJxx/vjZSUj2EfazlrC+88dLpR\niI7+C9u2rUNgYKCfovSvgwcPYubMedi37ygAoEqVchg0qC+qV6/u03JJolKlmjh2bDDI/nms/RsC\nA5/C+fPHi+STC/4khABJ17rru5q9XV0gW8Y5Gj16LPX6gR62jI/77Jncy5cvU6sNzLV8o7E3582b\n5/WyC6PY2FhGR99Ho7EKhRhFYDx1un7U6YLZqdNTBXrKxesuXbrEceMmsEGDtqxevREbNmzPDz+c\ncuPqx4ULF2g0hhLYkst5YaNON4CPP/6Mn2tzdzp48CCDgkpQiNm5tJCXU1HCGBsb6+9wiyTIDlyF\ny7lz56jXh9A+AIB7yVijGc4XX3zFJ/FZrVaGhZWlfcjGSwT+IRCfpfz86cldGKxatYoGQ3ECP2fz\nAXiVWu1gVqx4P69cueLvULNlsVj47LODqNMF0WDo7ahHLIEfqShPUa8P4qBBr3DMmLE0GJ5z4txM\noF4fzAsXLvi7agVSZmYmly1bxqZNOzAwsARNpjBWq/YgZ82a7ZXe3gcPHmSlSjUdXww/IrCVwDYC\n82g212N4eDlu2rTJCzWRsiOTcSHUoUM3CuHcZO/Z9YQ0GMJ56NAhn8SWkZHBJ554kkJE0j4fa1kC\nZgLNCHxPIabxvvsa+KTswiQpKclx/39THq3F/nz22YH+DvcOFouFDz3UkgZDVwL/5RD/BSpKO2q1\noY4P9rzPT72+Hz/8cLK/q1fgXLx4kffdV58mUx3a52Q+Q+ACgZU0mToyJCTK7XmOs7LZbNy4cSO7\ndOnFKlXqsWLFOmzduiuXLVvGzMxML9REyolMxoXQoUOHaDZH0Lnh67IuFipKU77wgm9axcnJyWzc\nuA2Nxodof1Qkw1FuGoHvCNxLjSaE+/bt80n5hcmsWbNpNHZy4nd2jnp9cIFrHT///Es0GLow76El\n0wk0JjDMyXN0sttjLxdVycnJrFq1DjWa17O5gnJ9+ZFms+++ZEu+504yliNw+VnlypWxYsUPMBp7\nwT7puzMSoCht0KJFCUybNskncT3xRG/s2BGK5OT1ADrjZicdLYAnAfwJtfohTJgwzSflFyZz5ixC\ncnJfJ9aMhEbTBMuXL/d5TM5KSEjAggVfIjV1BoC8Hv3RAFgAYC6AFCf2ngaDQetpiEXK/Plf4tSp\nEsjImAAgp/49XZCU9Cpee21MfoYm+ZlMxgVAo0aNsGnTaoSFDYXJ1AL2cWazGzLxGDSaYdDrK6NP\nnzpYuvQbnzw7uXfvXvz++1ZYLF8AOT6KrofF8h0WL16CEydOeD2GwsQ+N2wpp9ZNTy+N+Ph43wbk\ngvnzF0Clag8gwsktygGoC+DbPNc0m1egUaP67gdXxJDExImfISVlKHJOxNfXHYDVq1fhwoUL+ROc\n5HcyGRcQNWvWxJkzRzB79nO4//6PYDCUg9ncHXr9QBiNzyIwsDmMxnoYMMCGuLitmD59MlQq3/z6\npkyZgfT0AbC3gnNjgs32ND79dLZP4igszGYzAOcSrEbzn2P9gmHlyo1ISeno4laPA9iYxzq7ERBw\nHB07urrvouvixYu4cOEcgIedWDsIWm1TbNiwwddhSQWEx0PACCHaAvgI9mtcn5Oc4HFUdymdToce\nPXqgR48e2LdvH+Li4nDlyhXodDqEh4ejRYsWMBiyf97XmzZv3gGrtZ9T66ant0NsrG8ulRcWPXo8\niqNHv0Zqaqs81ryKjIxf0abN1HyJyxmJiUkAXH0WOBBC7HL0CcmuhXcWivIkJkwY65NRpgqrlJQU\nBASYkJbm3OOnNpsZqampPo5KKig8+ksRQqgBTAfQEsBZADuEEDEkD3gjuLtZ9erVfT5AQE4yMzOQ\nd6v4Og3S7eNj3rX69++LsWMrATgKoGKO66nVH6Nly9YFagCQkJAgONuqv+k/mM0JsFpbIzl5CIC2\nsH8XvwSVah50uo8xcuQQ9O/vzH30u0dYWBgyMi4DSAKQ93CnavVxREQ4e/tAKuw8vc5ZD8BRkido\nnxfuOwDeHQZKynflypUFEOfUukLEoXLlcj6Np6ALCwvD1KkToSitAGT3PdQGlepTBAfPwYwZH+Z3\neLl64om2MJkWubSN2bwIM2ZMwPTpPVGlytsICDBBpwuFVhuNxx8/iA0bYjBq1DAfRVx4mc1mNGnS\nEsA3Tqx9AEL8g+bNm/s6LKmgcLX7ddYF9ptHc7L83AvAJ7et48MO5JIv/PjjjzSbGzrx6IqNJlM1\nOYqPw+efz6OiFKPJ1JH2qTF/ohCTaDRWZpUqtfnPP//4O8Q7JCcnO0bU+sfJx5XiGBwcxfT09Bv7\nSEpK4qVLl255Tcre2rVraTSWz+V5bhKw0mDozBEjRvs7XMlN8MOjTfRwe6kAevTRR2EyXYAQX+e6\nnkr1KSIjdWjSpEk+RVaw9evXBxcvnsLUqY+iXbtlaNJkHnr1OoJff52HAwf+RPny5f0d4h0URcFr\nr73ieLQur8eVrkFRnsZbb70BjUZz41Wj0YiwsLBbXpOy16JFC/Tv/wSMxhYADmWzRjwMhqdx770J\nGDNmZH6HJ/mRRxNFCCHqA3ibZFvHzyMA2JilE5cQgmPG3HxerlmzZmjWrJnbZUr5Y9++fWjUqBUS\nE1+FzTYQt97jugq1eiqCg+dix44NiI6O9leYkhfYbDb07PkcfvllP5KT5wC4L5u1dsFo7IsePRpj\n9uxpEMK1MfClm0hi8uSPMXbsByBrIimpNYAAGAy7QcagW7cnMXPm1HzprCl5R2xsLGJjY2/8/M47\n74AuThThaTIOgP3rXQsA5wBsB9CDWTpwyVmbCq+jR49i4MDXsHnzRqhUjyI9PQxa7QVYrcvRsmVr\nzJjxIUqXLu3vMCUvIIlJk6Zi3LjJsFrLIzGxK+xzbV+G2bwYOt15jB79BgYPHiQTsZekpaXhhx9+\nwB9/7EJ6egaqVauAXr16IiwszN+hSR5yZ9Ymj6dQFEK0w81Hm+aSHHfb+zIZF3KnT5/GihUrcPXq\nVYSEhODRRx9FiRIl/B2W5AMZGRn4+eefsWpVLBISEhEaGohHHmmFDh06+GSAGUkqivySjPMsQCZj\nSZIk6S7iTjKWI3BJkiRJkp/JZCxJkiRJfiaTsSRJkiT5mUzGkiRJkuRnMhlLkiRJkp/JZCxJkiRJ\nfiaTsSRJkiT5mUzGkiRJkuRnMhlLkiRJkp/JZCxJkiRJfiaTsSRJkiT5mUzGkiRJkuRnMhlLkiRJ\nkp/JZCxJkiRJfiaTsSRJkiT5mUzGkiRJkuRnMhlLkiRJkp/JZCxJkiRJfiaTsSRJkiT5mUzGkiRJ\nkuRnMhlLkiRJkp/JZCxJkiRJfiaTsSRJkiT5mUzGkiRJkuRnMhlLkiRJkp/JZCxJkiRJfiaTsSRJ\nkiT5mUzGkiRJkuRnbidjIcTbQogzQojdjqWtNwOTJEmSpLuFJy1jAphCspZjWeWtoO5GsbGx/g6h\n0JDHyjnyODlPHivnyOPkO55ephZeiUKSJ7kL5LFyjjxOzpPHyjnyOPmOp8n4f0KIvUKIuUKIYK9E\nJEmSJEl3mVyTsRBijRAiLpulI4AZAKIB1ARwHsDkfIhXkiRJkoocQdLznQhRDsAvJO/L5j3PC5Ak\nSZKkQoSkS7dxA9wtSAgRSfK848cuAOK8EZAkSZIk3W3cTsYAJgghasLeq/o4gAHeCUmSJEmS7i5e\nuUwtSZIkSZL78mUELjlASO6EEG2FEAeFEEeEEG/4O56CSghxQgjxl+Mc2u7veAoSIcQ8IcS/Qoi4\nLK8Vc3TCPCyEWC2feMjxOMnPp9sIIUoLIX4XQvwthNgnhHjJ8bo8p26Ty7Fy6bzKl5axEGIMgESS\nU3xeWCEjhFADOASgJYCzAHYA6EHygF8DK4CEEMcB1CF52d+xFDRCiMYAkgAsuN6RUggxEcB/JCc6\nvuSFkBzuzzj9LYfjJD+fbiOEKAGgBMk9QggTgJ0AOgPoA3lO3SKXY9UNLpxX+Tk2tezIlb16AI6S\nPEEyA8B3ADr5OaaCTJ5H2SC5EUDCbS93BPCl4/9fwv4BcVfL4TgB8ry6BckLJPc4/p8E4ACAkpDn\n1B1yOVaAC+dVfiZjOUBI9koCOJ3l5zO4+YuUbkUAa4UQfwoh+vs7mEIgguS/jv//CyDCn8EUcPLz\nKQeOR1drAdgGeU7lKsux+sPxktPnldeSsRwgxG2yB53zGpKsBaAdgBcdlxwlJ9B+P0qea9mTn085\ncFx2XQLgZZKJWd+T59StHMfqB9iPVRJcPK88ebTpFiRbObOeEOJzAL94q9wi4CyA0ll+Lg1761i6\nzfXn2kleEkL8BPsl/o3+japA+1cIUYLkBSFEJICL/g6oICJ547jIz6ebhBAa2BPxVySXOl6W51Q2\nshyrhdePlavnVX71po7M8mOOA4Tcpf4EUEkIUU4IoQXwJIAYP8dU4AghFCGE2fF/I4DWkOdRXmIA\n9Hb8vzeApbmse9eSn093EkIIAHMB7Cf5UZa35Dl1m5yOlavnVX71pl4Ae1P9xgAhWe473PWEEO0A\nfARADWAuyXF+DqnAEUJEA/jJ8WMAgK/lcbpJCPEtgKYAwmC/lzcawM8AvgdQBsAJAN1IXvFXjAVB\nNsdpDIBmkJ9PtxBCNAKwAcBfuHkpegSA7ZDn1C1yOFYjAfSAC+eVHPRDkiRJkvwsP3tTS5IkSZKU\nDZmMJUmSJMnPZDKWJEmSJD+TyViSJEmS/EwmY0mSJEnyM5mMJUmSJMnPZDKWJEmSJD+TyViSJEmS\n/Oz/YRw6j5MarYIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105bac978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.scatter(x,y,s=(Cows+Sheep))\n",
    "plt.title(\"Size and Location of Farms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code in the cell block below generates one outbreak and returns a vector of the number of farms infected at each time point (up to 1000 time points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:61: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:94: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:103: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    }
   ],
   "source": [
    "#THE FOLLOWING FUNCTION RUNS THE OUTBREAK FOR 1000 TIMESTEPS (IDEALLY SO ALL FARMS ARE INFECTED)\n",
    "#BY THE 1000TH STEP. IT RETURNS A LIST OF THE NUMBER OF INFECTED FARMS AT EACH TIME STEP\n",
    "\n",
    "def Outbreaks(Size,N,Y0,farms,end,x,y,Cows,Sheep,Maxtime=1000):\n",
    "    #This is an attempt of converting the Matlab Program 7.6 Code into Python\n",
    "    import numpy as np\n",
    "    import pandas as pd\n",
    "    from math import pi\n",
    "    #Generate dataset (assume there are the same amount of animals at each farm)\n",
    "    #np.random.seed(10)\n",
    "    #x = Size*np.random.rand(N)\n",
    "    #np.random.seed(11)\n",
    "    #y = Size*np.random.rand(N)\n",
    "    #Cows = np.array([250]*N)  #Number of cows at each farm\n",
    "    #Sheep = np.array([250]*N)  #Number of cows at each farm\n",
    "    Status = np.array([0]*N)    #Initial Status of each farm\n",
    "    init_ind = np.random.randint(0,N)\n",
    "    for i in range(Y0):\n",
    "        Status[init_ind] = 6\n",
    "\n",
    "    #Cows are 10.5 times more susceptible to disease than sheep\n",
    "    Suscept = Sheep+10.5*Cows\n",
    "    Transmiss = 5.1e-7*Sheep + 7.7e-7*Cows \n",
    "\n",
    "    #Calculates which grid square a particular location is in (turn a 2-d coordinate into a scalar)\n",
    "    def WhichGrid(x,y,XRange,YRange,XNum,YNum):\n",
    "        #Essentially: floor(Unif[0,1)griddim)griddim+floor(Unif[0,1)griddim)+1\n",
    "        #Returns a number from 1 to griddim^2\n",
    "        return(np.floor(x*(XNum/XRange))*YNum+np.floor(y*(YNum/YRange))+1)\n",
    "\n",
    "    #Set up the grid\n",
    "    grid = WhichGrid(x,y,Size,Size,10.0,10.0)\n",
    "    tmp = sorted(grid) #Sort grid values\n",
    "    i = np.argsort(grid) #get indexed values after sort\n",
    "    x = x[i]\n",
    "    y = y[i]\n",
    "    Status = Status[i]\n",
    "    grid = grid[i]\n",
    "    Xgrid = []\n",
    "    Ygrid = []\n",
    "    Num = []\n",
    "    first_in_grid = []\n",
    "    last_in_grid = []\n",
    "    Max_Sus_grid = []\n",
    "    index_inf = np.where(Status==6)[0].astype(int)\n",
    "    for i in range(1,int(max(grid))+1):\n",
    "        #turn the grid square number into an x-coordinate and y-coordinate (should not exceed XNum)\n",
    "        Xgrid.append(np.floor((i-1)/10))\n",
    "        Ygrid.append((i-1)%10)\n",
    "        m = np.array(np.where(grid==i))\n",
    "        Num.append(m.shape[1])\n",
    "    \n",
    "        if Num[i-1] > 0:\n",
    "            first_in_grid.append(m.min()+1) #Add the \"+1\" here so the indicies match those in the Keeling code\n",
    "            last_in_grid.append(m.max()+1)\n",
    "            Max_Sus_grid.append(Suscept[m].max())\n",
    "        else:\n",
    "            first_in_grid.append(0)\n",
    "            last_in_grid.append(-1)\n",
    "            Max_Sus_grid.append(0)\n",
    "\n",
    "    #Work out grid to maximum grid transmission probabilities\n",
    "    from numpy import ndarray\n",
    "    MaxRate = ndarray((max(grid),max(grid)))\n",
    "    def Kernel(dist_squared):\n",
    "        dist_squared = np.asarray(dist_squared)\n",
    "        is_scalar = False if dist_squared.ndim > 0 else True\n",
    "        dist_squared.shape = (1,)*(1-dist_squared.ndim) + dist_squared.shape\n",
    "        K = 1 / (pi * (1 + dist_squared**2))\n",
    "        K[(dist_squared < 0.0138)] = 0.3093\n",
    "        K[(dist_squared > 60*60)] = 0\n",
    "        return(K if not is_scalar else K[0])\n",
    "\n",
    "    #Determine maximum number of animals to be infected in each grid square\n",
    "\n",
    "    for i in range (1,int(max(grid))+1):\n",
    "       for j in range(1,int(max(grid))+1):\n",
    "           if ((i-1)==(j-1)) | (Num[i-1]==0) | (Num[j-1] == 0):\n",
    "               MaxRate[i-1,j-1] = np.inf\n",
    "           else:\n",
    "               Dist2 = (Size*max([0,(abs(Xgrid[i-1]-Xgrid[j-1])-1)])/10)**2+(Size*max([0,(abs(Ygrid[i-1]-Ygrid[j-1])-1)])/10)**2\n",
    "               MaxRate[i-1,j-1] = Max_Sus_grid[j-1]*Kernel(Dist2)\n",
    "\n",
    "    #Susceptible, Exposed, Infectious, Reported.==> latent period is 4 days\n",
    "    i=1; S=len(np.where(Status==0)); E=len(np.where(np.logical_and(Status>0, Status<=5)));I=len(np.where(np.logical_and(Status>5, Status<=9))); R=len(np.where(Status==10)); R2=len(np.where(Status>9)); CullSheep=0; CullCattle=0;\n",
    "    i=i+1;  IterateFlag=1;\n",
    "\n",
    "    def Iterate(Status, x, y, Suscept, Transmiss, grid, first_in_grid, last_in_grid, Num, MaxRate):\n",
    "        Event = 0*Status\n",
    "        INF = np.where(Status>5)[0]\n",
    "        NI = INF.size # Note reported farms still infectious\n",
    "        IGrids = grid[INF]-1\n",
    "        \n",
    "        for ii in range(NI):\n",
    "            INFi = INF[ii]\n",
    "            trans = np.multiply(-Transmiss[INFi],Num) #transmissibility of infected farm to all other grid squares \n",
    "            maxr = MaxRate[IGrids[ii],:] #max number of animals to be infected in infected grid square\n",
    "            # Elementwise multiplication\n",
    "            rate = np.multiply(trans, maxr) #max number of animals to be infected in each grid square based on infected grid square\n",
    "            MaxProb = 1 - np.exp(rate) #Max probability that infected farm infected noninfected farm\n",
    "            rng = np.random.rand(len(MaxProb))\n",
    "            m = np.where((MaxProb - rng)>0)[0]  #these grid squares need further consideration\n",
    "            for n in range(len(m)):\n",
    "                s = 1\n",
    "                M = m[n]\n",
    "                PAB = 1 - np.exp(-Transmiss[INFi]*MaxRate[IGrids[ii],M]) #Max probability that infected farm infects noninfected farms under consideration\n",
    "                if (PAB == 1):\n",
    "                    # Calculate the infection probability for each farm in the susceptible grid\n",
    "                    leng = last_in_grid[M]-first_in_grid[M]+1\n",
    "                    R = np.random.rand(leng)\n",
    "                    for j in range(leng):\n",
    "                        ind1 = first_in_grid[M]+j-1\n",
    "                        Q = 1 - np.exp(-Transmiss[INFi]*Suscept[ind1]*Kernel((x[INFi]-x[ind1]**2+(y[INFi]-y[ind1])**2)))\n",
    "                        if ((R[j] < Q) & (Status[ind1] == 0)):\n",
    "                            Event[ind1] = 1\n",
    "                else:\n",
    "                    R = np.random.rand(Num[M])\n",
    "                    # Loop through all susceptible farms in the grids where an infection event occurred.  \n",
    "                    for j in range(Num[M]):\n",
    "                        P = 1 - s*(1 - PAB)**(Num[M] - j)\n",
    "                        if (R[j] < (PAB / P)):\n",
    "                            s = 0\n",
    "                            ind1=first_in_grid[M]+j-1\n",
    "                            Q=1-np.exp(-Transmiss[INFi]*Suscept[ind1]*Kernel((x[INFi]-x[ind1]**2+(y[INFi]-y[ind1])**2)))\n",
    "                            if ((R[j]< Q/P) & (Status[ind1] == 0)):\n",
    "                                Event[ind1] = 1\n",
    "        # Evolve the infection process of those farms which have been exposed and already infectious ones.  \n",
    "        Status[Status > 0] += 1\n",
    "        Status = Status + Event\n",
    "        #m=np.where(Status==13); # Initiate Ring Culling Around Reported Farm\n",
    "        #for i in range(len(m)):\n",
    "        #    Status[m[i]]=-1;\n",
    "        return {'Status':Status,'NI':NI}\n",
    "    \n",
    "    S=[]\n",
    "    E=[]\n",
    "    I=[]\n",
    "    R=[]\n",
    "    R2=[]\n",
    "    CullSheep=[]\n",
    "    CullCattle=[]\n",
    "    t=[]\n",
    "    t.append(0)\n",
    "    results = np.c_[np.array([1]*N),np.arange(1,N+1),np.array([0]*N)]\n",
    "    Maxtime = 0\n",
    "    while(Maxtime<1000):\n",
    "        Status=Iterate(Status, x, y, Suscept, Transmiss, grid, first_in_grid, last_in_grid, Num, MaxRate)['Status']\n",
    "        Sus=np.where(Status==0)[0]; Exp=np.where(np.logical_and(Status>0, Status<=5))[0]; Inf=np.where(Status>5)[0]; \n",
    "    \n",
    "        S.append(len(Sus)); E.append(len(Exp)); I.append(len(Inf)); \n",
    "        t.append(t[i-2]+1);i+=1;\n",
    "    \n",
    "        Maxtime+= 1\n",
    "    vec = np.c_[list(range(1,len(I)+1)),I]\n",
    "    return(vec[:,1])\n",
    "\n",
    "test = Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=20,x=x,y=y,Cows=Cows,Sheep=Sheep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the above algorithm is implemented 100 times, so the number of infected farms can be averaged over 100 simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:61: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:94: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entering for Loop\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:103: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from random import randint\n",
    "Num_outbreaks = 99\n",
    "num_inf = Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=20,x=x,y=y,Cows=Cows,Sheep=Sheep)\n",
    "print(\"Entering for Loop\")\n",
    "for i in range(Num_outbreaks):\n",
    "    num_inf += Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=20,x=x,y=y,Cows=Cows,Sheep=Sheep)\n",
    "    print(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot our results. Two graphs are presented: the first one is a trajectory across all 1000 timesteps, the second graph takes a closer look at the beginning of the outbreak (recall that we want to implement an intervention, so ideally an action would occur towards to beginning of the outbreak). Based on the second graph, it seems that it would be interesting for data to be collected on day 20 (an average of 3 farms out of the 75 are infected) and an intervention can be implemented on day 35 (an average of 8 farms are infected). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1090cdc88>"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAADrCAYAAAD5e5gnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVNX5x/HPlyqgNOWHHdSIYMGoWGIsE3tFjT0JYmLU\nRJGosWBJQE000SR2Y4xRsaOIimiMSFhjVxQVJaBGLKCgCAKC9Of3xznD3h1ndmeXKbs7z/v1uq+9\nc+szd2bvM+fcc++RmeGcc85VkhblDsA555wrNU9+zjnnKo4nP+eccxXHk59zzrmK48nPOedcxfHk\n55xzruJ48isSST0lrZRUkmMs6UNJexdhuylJnxRoWyslbRrH75B0WSG2W679NDaSfifpC0mfljmO\nnN/FSvo8ACT9VdLF5Y7DfZsnv0jSiZImSVoo6TNJN0nqVI/1P5S0VzFjrIPFod6SyaKEGhwvgKRN\nYtw3NXQ/8TNfIWlBYriuoTEVmqRhku7Kc9mNgbOB3ma2/mrssxA/dmr7bFfrc0+SdFn8n10maWiW\n+T+S9JGkryU9LKlLYl5bSbdJmhf/38/KWPe7kl6L54MJkrbNEcM/E9+dpZKWJF7fZGa/NLPfFeL9\n1kXShZI+iPv+RNL9iXlVkk4qRRxNhSc/QNKvgT8AvwY6ArsAPYCxklrnuRkDtBoxtGzouk1Yg48X\ncALwNnCspDarsZ/nzWytxDC4PkFIalWf5eupPkliY+BLM/uyWMEU0Op87knvAecCj5NxrCRtBdwM\n/BjoDiwCkj+UhgGbEY7bD4DzJO0f120DPArcCXQGhgOPZjsXmNmB6e8OcA/wx8R36bQCvc86SRoI\n/ATYO8bSD3g6GWqpYmkqKj75SepI+EcYZGZPmdkKM/sIOAboSfhCfau6JvkLOf463xh4LP7qOiex\ni5MkzZD0aUyy6fWHSRop6S5J84CBkjpJ+kdcdnr8ZdsiLr+ZpH9Lmh2rtu7OVTKV1Cf+Ajy2Acej\nXXyvcyS9A+yYMX99SQ9J+jzu44zEvJ0kvShpbnwP19fx48Hiem9LOiSxndbxfeb6tS1gAOFz+xI4\ntL7vM7m5LNv/qaTJkuZL+p+kUxLzUvGzOU/SZ8BtkoZKejB+lvMlvSVpc0kXSJoVSx/7JrZxYtzu\n/HgMf5RPbAol3VMlvRuP8Q1x+j7AU8D68ft3W5y+i6QX4rJvSNozsa2ukm6P3805kkZJag/8M7Gd\n+ZLWVTBE0vvxcxmhmqWoAfE9zpZ0YR7HfB1JT8XtVymUWpF0o6Q/Zbzn0ZLOzLYRM7vTzJ4EFmQe\nK0LSG21mz5nZQuA3wA8ldYjzTwAuM7N5ZjYFuAU4Mc5LAS3N7FozW2Zm18ft51Ozk/mZrTpvJL47\n58b/n08lHS7poPiZfilpSGLdWo97hn7Av8xsWjw2s8zs1rid3wO7AzcoUbshqbeksXG/UyQdnRH3\nzdk+pzj/6vjdnhe/71vlcWwalYpPfsCuwBrAqOTE+A/zBJA+aeWsrjGzAcDHwCHxF1/yHzgFfAfY\nDzhfNa+F9AceNLNOwL3AHcBSwi/S7eI6P08s/3tgPaAPsBHh5F+DpO2BJwnJfEScdqOkG2s5BklD\ngU2ATYH9gYFUJ6kWwGPARGB9YG/gTEn7xXWXA78C1ga+F+fn8+t3OPFHRnQQMMPM3syx/G6EX/NP\nAA/GGAtpFnCwmXUEfgpcLWm7xPzuQBfCD55TCCe8QwglhS6E4zM2Lrs+cBnwN4B48r0WOCBu/3vA\nG/WI7WDCia4vcIyk/c3saeBA4NP4/fuZpA2AMcClZtYFOAd4SNLacTt3Eb73WwL/B1xtZouAAxLb\n6WhmM4HBhO/qHoTv31zgxvh+tiSUqH4c3+vawIa1xK+47KXAOvG93xPn3QEcH3/cIGkdwnfonm9v\npk5bAqu+P2b2AbAE6BUTyHrJ+cBbQPoEvlV8nfRmYn59ZJ43ugNt4/5/C9xKOB7bERLUbyX1iMvm\nPO5ZvAScIOkcSf2UqEkys4uAZ4HT07Ub8Xs4Frgb6AYcB9wkqU9imz8iy+ekUELeHdg8nruOJvwI\nbVI8+YUPdraZrcwybybhnzmtIdU1l5jZN2b2NnA7cHxi3gtmNjqOdyKcwM6Ky38BXEP4UmJm/zOz\ncfGX6GzgamBPatqTUF0zwMyeSE80s9PN7PQ84z0a+L2ZfWVm0wkn6vT73hFYx8x+Z2bL46/MWxMx\nvm5mr5jZylh6viVLjEnp7d4DHCxpzfh6AOHknMtA4DEzW0xIfgdI6pbn+8u0SywZpYedzOyJxC/o\n/xBKVbsn1lkJDI2fxeI47T9mNtbMVgAjCd+bP8TXI4CeCrUM6fW3kdQu/kKfXI94/2Bm883sE2A8\n8N04PfO7+RPgiVgyIibICYTjvB4hyf0ilnyWm9mzObYDcCpwsZl9ambLgEuAo+IJ9ijCZ/GcmS0l\nlLCy/S8ljUksfxHwPUkbmNmrwDxCwoPwvRof/xfqa824raT5wFpxHhnz0/PqWrchksd0GeH/K/29\n6ApcY2YL4/dgMpCu8fgF2Y/7t87bZnYPcAbhB2sVMEvSebXEcQgwzcyGx//XNwgFgKMTy2T9nAg/\n0NcC+khqYWZT44+kJsWTH8wmVMNkOxbrxfmrI9l44GPCr+O06YnxHkBr4LP0iZhwzaIbgKTuku6P\n1SbzCMkhMzGfSriG9Z/ViHf9LDEnY1w/mSyACwglByT1kjRGoQHBPEJJNRljVmb2KfA84R+7M+HE\nnPXXvqR2hBPug3HdN4APCb9SG+IlM+uSGF6RdKCkl2J10FxCSTT5Pr6IJ4SkzxPj3xB+UFniNcCa\nsUbhWMKJ7dN4vLaoR7zJk8wiqk/kmXoAR2d8Vt8H1iXUGswxs8wTfC49gYcT25lMKOV3J/yPrPoe\nx9JjbaUAy1h+ITCH6v+LO6muBfgJtf8Iqs3XhB+USZ0IVaRfx9cds8xLr5ucl54/v4GxJH2Z5Xsx\nKzH/G6o/0x7kPu7fYmb3mtm+MdZfAJcpUd1OzRJoD2DnjO/HjxLbzvk5mdl44AZCKXSWpL9JaugP\ng7Lx5AcvEqpDjkxOjKWQA4BxcdJCoH1ikXUztpPrgvLGGeMzcqzzSYxj7cSJuJOZbRPnXw6sALaO\nVQ0DqPn5GSH59ZD0lxyx5OOzLDEnY5yWkSw6mln6et1fCf+g34kxXkT+37F01efRhBLxZzmWO4Jw\nYvpbTLKfEU7mBan6lNQWeAi4Evi/WGX4BDV/NWd+1vVqTGDh2vJ+hO/QFODvuRbNZ3M5pn8M3JXx\nWa1lZlcSPseuyn7NONv2PiZU0ya31T7+aEkffwAUrhvW9YMnufyahNJP+vaMu4HDFK739gYeqWNb\nueJ+h+oSFJI2A9oA75rZ3Bj3dxPLb0toQJVet2/G9vrG6Q3R0MYmuY57rv+NsLPQbmEkoep26xwx\nfAw8k+X7ka4hErV8TmZ2vZn1I1Qv9yI0PGpSKj75xV+/lwDXS9pfobFFT+ABwkki/cvzDeAgSV0k\nrQtkXoSfRbhWl+lihUYkWxEuqI/IEcdnhOq1v0haS1ILhUYue8RF1iQk4Pmx6iHbl20BIWHvIemK\nPN5+Ng8AF0jqLGlDQlVK2ivAAoXGHu0ktZS0taR+iRgXAIsk9QZ+Wct+MqvXHga2J1znuLOW9QYC\n/yD8U28bh+8D20raOsvy9a2qbhOH2cBKSQcSrr3WJu99SPo/SYfFay7LCJ/pigZuV7UsczdwqKT9\n4ue0hkKDiw3id+2fhGs8neN3Pv09mwWsnaiihVADcbmqG6Z0k9Q/zhsJHCLp+wqtJC+l9vOKCP9H\n6eUvA140sxkAsap9AuE7MNLMluTckNRK0hpAS6B1fI/pfd8T3/9u8VhfBjwUSzDE7V8c338fwrX1\nO+K8KmCFpMEKt0QMJlTl/ruW95V+b9mmNbR1a23HveZOpIEKDWfS544DCdcoX46LZJ6fxhCuf/4k\nfv6tJe0Y/2/Tsn5OCtcUd1ZozLYIWEzu73CjVfHJD8DMrgIuBP5EqOt/CfiI0Gx4WVzsLsJF7w8J\nDUrup+avqSsI/0xzJZ1N9YXuZ4D3Cc2Or4rXXkjMTzqBcOKdTKhieJDqEuYlhOQwj9Do5KEs66eT\n+b7AgZIuAVBotfXX2g5BYvyS+N6nxfd5Z3p+vE5xCOEX8wfAF4TreukT5TmEqpP5cXrmMcocX/U6\nXjsbRahiq9H4KC0m/b0I10g+Twyvx1hPyPHe8r7nzMwWEBLwA4TP4HjCddTM9eraR67XLYCzCDUA\nXxKuJeb6kZC53br2mzye04HDCN/rzwm/9H9N9f/8AELynUI4MQ6O600B7gM+UGgFui7huu9o4ClJ\n8wm1JTvF5ScDpxMabH1KOGa13SdohMQ0NL7/7ajZ2AlCLcA21F3leSvh5HscoZZhUXpbMa5fxH3N\nAtpRs/HVUOB/hO/6eMItCk/FdZcChxO+T3Pj38PNbHkd8eT6HtT1GeaS87hnMZ/wWX8UY/4D4Zru\nC4ltHRU/02vM7GvCj7rjCN/FzwjnsPRtQ0b4TLN9Th0J/99zCOfD2cBVtbyPRklWxM5sJV1AOGAr\ngUmElnMdCKWfHoQDd4yZfVW0IFyTIek3hBZk2ZKYqxCSdgfuNrMedS7sikLS7cB0M/tNuWMplqKV\n/GLV4cnA9vG6VUvCr4whwFgz60W4njYk1zZc5ZDUFfgZ4Relq1CxKu1Mcl8HdaVRqAcRNFrFrPac\nT6hWaa/wFIz2hGqR/oRqDeLfw4sYg2sCJJ1MqJb7p5k9V+54XHnEa29zCS0OrylzOJWuYI+ha6yK\nXe15CvBnQvPdf5nZAElzYwu69JM65qRfO+ecc6VQtOcSxqbFZxIaMMwDHpRU48K2mZmkXA8dbta/\nOpxzztWfmRWkSraY1Z79CPdrfRlbSY0iPMppZmxBhsKTJj7PtQEz86Gew9ChQ8seQ1Mc/Lj5cfPj\n1viHQipm8ptCeHRUu1i9uQ+hCf9jVN+QPJD8b2J1zjnnCqJo1Z5m9qakOwk3rK4EXie05FsLeECh\nb6kPCb0nOOeccyVTzL7IsPAopSszJs8hlAJdEaRSqXKH0CT5cWsYP24N48et/Ira2nN1SLLGGptz\nzrnSk4Q1gQYvzjnnXKPkyc8551zF8eTnnHOu4njyc845V3E8+TnnnKs4nvycc85VHE9+zjnnKk69\nkp+klpI61r2kc845VxhmsGJFYbdZ5xNeJN0HnAqsAF4FOkm6Nj69xTnnXIGYwcqV4US/YgUsX149\nnny9ZAksXhz+Ll0apmcO6eVzDZnzV6yo3veyZdXbXrwYFi4M09JxJGNMD0uXhnWS82obT79Oj9c2\nALQocD1lnU94kfSmmW0r6cfA9oSe11+30Dt70fgTXpxzkP1EmxyWLIFFi8LJd9myMKTHFy2CBQtq\nnrhzDekT/TffhGHx4uwn4nSCqmt65gk/WyJbuhTmzQtxLlsWpknQsiW0ahX+pofk6zXWgLZtw982\nbaB16+pl0kN62datq18nx5PLJZeXwjbbtg1/11gD2rcP48l4WrSo+Tq9Tua8XOPp95me1qJF7kHx\nmS6FfMJLPs/2bCWpNaHH9RvNbJn3tedc05JZosgckvPSSSA95CpZLF9enXiSwzffhHWSQ/rEvnx5\ndYLKHNKJJ73fJUuqq7qSJ9nMoW3b6pNz69Y1hw4dYK21qpNDbUObNmE7XbtCu3Zhu61a1X1SzjU9\neULPTDDp8TZtoFOnsN90Yip0Ccdll0/y+xuh94W3gP9I6knonNY5lwezcDJftChUH2X+/eqrMD8z\nGSxdGpJBOqmk11m0KCyfq+oqXaJIJpLFi6tPyNmG5C/ydIkiPaQTSbJUkR7SCaZDh3AC79QJuncP\n20iXHFq3Dn+TpY9sQ7t21fts1666pOHJwBVDvR9sHfvma2mhg9q6lt0CuD8xaVPgN8DdwAigB7Fb\nIzP7KmNdr/Z0JWNWfc0inSzmzoUvvoAvvwyJJ9uQTkrJ1+lqs2++qZ6WLlWkk0QyYXTuHE72ySST\nHk+ukx46dKgulWRWXaXHkwksXYrxJOKaukJWe+Zzza8LcALQk+qSopnZ4HrtSGoBzAB2As4AZpvZ\nlZLOB7qY2ZCM5T35udW2YkVIXpMmwf/+Bx99BHPmhOtAM2fCrFnVCU6qeS2lSxfo1i1Ug6255rcT\nV7ahffuQbNJDhw7hb6uidh7mXGUo9TW/J4AXCdWeKwEBDclK+wDvm9knkvoDe8bpw4EqQkMa5/Ky\nYAHMmAHTp8Mnn4TxmTNDae3zz2HaNPj001CS69wZttoKNt8cevSAbbYJyWzddUMVXbdusM46oXTm\nnKsM+SS/tmZ2dgH2dRxwXxzvbmaz4vgsoHsBtu+agZUrQ/KaMaN6mD79269XrIANNgjDRhuFv1ts\nUV1a69kzTOvQobqlmHPOpeWT/O6VdArwGLAkPdHM5uS7E0ltgEOB8zPnmZnlaj06bNiwVeOpVMp7\nP24m5syB//4X3nwT3n+/uuSWLr116lSd2DbcMPzdY4+arzt18qTmXHNXVVVFVVVVUbadzzW/QcDv\nga8I1Z4Qctamee9EOgz4pZkdEF9PAVJmNlPSesB4M+udsY5f82vCzML1tMmTQ6JL/l20CPr0gb59\nq6si08lu/fXDNTfnnMtU6mt+vwY2M7PZq7Gf46mu8gQYDQwE/hj/PrIa23ZlZgbvvBOS29ix1eMt\nWsCWW4ahTx84/PDwd4MNvNTmnCuvfEp+TwFHmNnCBu1A6gB8BGxiZgvitK7AA8DG+K0OTY5ZKMGN\nGwdVVfDSS6F15JZbwv77w7bbhiTXrVu5I3XONSelvtXhEWArYDzV1/zqfatDvQPz5NeoTJ8OTz8d\nhnHjQvP9ffaB3XaD3XcPDUy8NOecK6ZSJ7+B1Ly9QYTkN7wQAdSyX09+ZZQu3Y0ZA6NHw9SpsPfe\nIeHtvTdsmvcVX+ecK4ySJT9JrYCnzSxViJ3Vhye/0lu4EJ55Bp58MiS9FSvgkEPCsPfefh+cc668\nStbgxcyWS1ohqXPmNTnX9JnB88/D44+H63avvQbbbw/77guPPBJuBveqTOdcc5RPtedoYDtgLJBu\n9OLX/Jqwjz6CO+6A4cPDtbsjj4Rdd4VddglPQ3HOucao1Lc6jIpDkmelJmjSJLjgglDKO/54ePDB\nUNLz0p1zrtLUu1eHUvGSX+HMmAG//W24jnfRRXDqqX4juXOu6Slkya/OTk4k9ZI0UtJkSdPi8EEh\ndu6Ka968UNLr2zfcczd1Kgwe7InPOefyqfa8HRgK/AVIAT8FWhYxJlcAU6bA0UeHas2JE2Hjjcsd\nkXPONR75dG/ZzsyeJlSRfmRmw4CDixuWWx0PPhhuPB80KDRs8cTnnHM15VPyWyypJfB+fMj1p0CH\n4oblGuLrr+HCC+Ghh8JTWPr2LXdEzjnXOOVT8jsTaA8MBvoBPyE8jNo1IhMmQL9+MH9+aNXpic85\n53LL2dpT0lNmtl8cv9DMLi9pYN7aMy8rV8I558C998LVV4dbGJxzrjkq1X1+yWfyHw2UNPm5upnB\n2WfD66/D22/DOuuUOyLnnGsa8rnm5xqhlSvhjDPCDevjxvmTWZxzrj5qS36bxkebCdhE0mOJeWZm\n/fPZgaTOwK2EbpGMcKvEe8AIoAc5+vNzuU2bBiecAK1awfjx0LFjuSNyzrmmpbZrfqla1jMzeyav\nHUjDgWfM7LbYS0QH4CJgtpldKel8oIuZDclYz6/5ZfHuu6GHhcGD4de/Dr2lO+dcJShpf36rtXGp\nEzDRzDbNmD4F2NPMZklaF6gys94Zy3jyy5DuU+/SS+FnPyt3NM45V1olfbzZatoE+ELS7ZJel/R3\nSR2A7mY2Ky4zC+he5DiavJdfDonvd7/zxOecc6ur2A1eWgHbA4PM7FVJ1wA1qjfNzCRlLeINGzZs\n1XgqlSKVShUv0kbspZegf3+45RY4/PByR+Occ6VRVVVFVVVVUbZd7GrPdYEXzWyT+Ho34AJgU+AH\nZjZT0nrAeK/2zO6tt0LnsrffDgcdVO5onHOufEpyn19m605Cq89Vr/Np7RmT2yeSepnZu8A+wDtx\nGAj8Mf59pCHBN3dTp8IBB8D113vic865QsqntecRwLrA3YQEeDwwy8zOzGsH0raEWx3aAP+juleI\nB4CNyXGrQ6WX/KZPD72rX3opnHhiuaNxzrnyK2lrT0mvmdkOdU0rtEpOfvPmQSoVHlV23nnljsY5\n5xqHUrf2bC9ps8TONyU86NoVwcKFoYpzt93g3HPLHY1zzjVP+ZT8DgBuAabFST2BU8zsX0UNrAJL\nfmahA9oOHUIDF7+B3TnnqpXqwdYAmNmTknoBW8RJU8xsSSF27mq64Ybw6LIXXvDE55xzxZRPya8D\ncDawsZmdLGlzYAszG1PUwCqs5Pfqq3DwwfDii7DZZnUv75xzlabU1/xuB5YCu8bXnwK/L8TOXTB3\nLhx7LPz1r574nHOuFPJJfpuZ2R8JCRAzW1jckCqLWXhc2SGHwJFHljsa55yrDPk83myJpHbpF7Hl\np1/zK5Brr4UZM+D++8sdiXPOVY58kt8w4ElgQ0n3At8HTixiTBXjlVfg8svDQ6vbti13NM45Vzny\neranpHWAXeLLl83si6JGRfNv8PLNN7DttiH5HXVUuaNxzrnGr9RPeBlnZnvXNa3QmnvyO/VUWLAA\n7r233JE451zTUKoHW7cjPMmlm6SuiVkdgQ0KsfNKNWIEVFWF2xucc86VXm3X/E4FfgWsD7yWmL4A\nuKGYQTVn8+aF53XedRd07FjuaJxzrjLlU+15hpldX6J4kvttdtWe6ceXde8ON95Y7micc65pKfVN\n7iapS2LnXSSdlu8OJH0o6S1JEyW9Eqd1lTRW0ruSnpLUuQGxNzk33hgeX/bnP5c7Euecq2z5lPze\nNLNtM6a9YWbfzWsH0jRgBzObk5h2JTDbzK6UdD7QxcyGZKzXrEp+r70GBx7ojy9zzrmGKnXJr4Wk\nVctJagm0rud+MoPtDwyP48OBw+u5vSZlyRIYMACuu84Tn3PONQb5JL9/AfdL2lvSPsD9hJve82XA\n05ImSDo5TutuZrPi+Cygez221+RcdRV85zvh+Z3OOefKL58nvJwPnAL8Mr4eC9xaj31838w+k9QN\nGCtpSnKmmZmkrPWbw4YNWzWeSqVIpVL12G3j8N57cM018PrroIIU1p1zrjJUVVVRVVVVlG3n+4SX\n9oQujabUuXDt2xkKfA2cDKTMbKak9YDxZtY7Y9kmf83PDPbdN3RVdNZZ5Y7GOeeatpJe85PUH5hI\nrOqUtJ2k0flsXFJ7SWvF8Q7AfsAkYDQwMC42EHik/qE3fnffDXPmwBlnlDsS55xzSfm09nwd2ItQ\nOtsuTnvbzLauc+PSJsDD8WUr4B4zuyI+MeYBYGPgQ+AYM/sqY90mXfL78kvYaisYMwb69St3NM45\n1/SV5PFmCcvM7CvVvGC1Mp+Nm9k04Fu3RMTbHvbJK8Im6rzzQgMXT3zOOdf45JP83pH0Y6CVpM2B\nwcALxQ2raXv2WXjqKXjnnXJH4pxzLpt8bnUYBGxF6MD2PmA+cGYxg2rKli4NPTZcc40/u9M55xqr\nnNf8JN1lZgMknWlm15Q4riZ7ze+KK+D55+Gxx/zWBuecK6SS9OcnaTLhutyTQCpzfvJxZcXQFJPf\nBx/ATjvBhAnQs2e5o3HOuealVA1ebgbGAZtSs0sjCE9t2bQQATQXZjBoEJx7ric+55xr7PK51eFm\nM/tFieJJ7rdJlfxGjoRhw2DiRGhd3yefOuecq1NJqj0zdtiS8PzNVSVFM/u4EAHUss8mk/yWLIHN\nN4d77oHddy93NM451zyV9D4/SWcAQ4HPgRWJWdsUIoDm4LbbYJttPPE551xTkU+15/+Anczsy9KE\ntGq/TaLkt3Rp6LFh5MjQ2MU551xxlLo/v48J9/a5LG69Fbbc0hOfc841JfmU/G4DegGPA0vjZDOz\nvxQ1sCZQ8vvkE9h+exg3Dvr2LXc0zjnXvJX62Z4fx6FNHES41aHiXXwxnHaaJz7nnGtq8mrtWQ6N\nveQ3Y0Zo5DJ1KnTrVu5onHOu+StJyU/SY7WsZ2bWP58dxNskJgDTzezQ2J3RCKAHObozagpuvhkG\nDPDE55xzTVFtjzdL1bKemdkzee1AOhvYAVjLzPpLuhKYbWZXSjof6GJmQ7Ks12hLfrNnh0YuzzwD\nffqUOxrnnKsMJb/JvcEblzYE7gB+D5wdS35TgD3NbJakdYEqM+udZd1Gm/wGDQoPrb7++nJH4pxz\nlaPUDV5Wx9XAuUCyc5/uZjYrjs8iPDmmyXjvPbj/fpgypdyROOeca6iiJT9JhwCfm9nEXFWoZmaS\nchbvhg0btmo8lUqRSmXdTEldfDGcfTass065I3HOueatqqqKqqqqomy7aP35SbocGAAsB9YglP5G\nATsCKTObKWk9YHxTqfacNAn23x/efx/aty93NM45V1lK9YSXHSStD/xMUtfMoa4Nm9mFZraRmW0C\nHAf828wGAKOBgXGxgcAjq/smSuXPf4YzzvDE55xzTV0p+/NLF+P+ADwg6STirQ713E5ZTJ0Kjz8O\nV19d7kicc86tLu/PL08DBoTbGy64oNyROOdcZSpHf37bAnsQSm/Pmtmbhdh5HftsNMnv449h223D\ntb611y53NM45V5lK2quDpF8B9wDdCLcl3C1pcCF23lScdhqcc44nPuecay7yqfacBOxiZgvj6w7A\nS2ZW1M5sG0vJ77XX4LDD4IMPoE2bckfjnHOVq9T9+QGszDHe7F11FZx1lic+55xrTvK5yf124GVJ\nowjdGR0O3FbUqBqJDz6Ap5+GW24pdyTOOecKKd8GLzsAu1Hd4GVi0QNrBNWep58OnTrB5ZeXNQzn\nnHM0oQdbr45yJ7/Zs6FXL5g8GdZdt2xhOOeci8pxza/i3HorHH64Jz7nnGuOvOSXxYoVsOmmMGoU\n7LBDWUJwzjmXoWQlP0mtJI0vxI6akjFjYP31PfE551xzVWvyM7PlwEpJnUsUT6Nw442hsYtzzrnm\nKZ9bHRYI2jJlAAATR0lEQVQCkySNjeMQuuJrlk95mToV3nwTjj663JE455wrlnyS36g4pC/AKTHe\n7Nx0E/z859C2bbkjcc45Vyz53ufXHtjYzKbkvWFpDeAZoC3QBnjUzC6IfQGOAHoQuzQys6+yrF/y\nBi9ffw09esDEibDxxiXdtXPOuTqU+sHW/YGJwJPx9XaSRte1npktBn5gZt8F+gI/kLQbMAQYa2a9\nCP0FDlmN+Avqnntgjz088TnnXHOXz31+w4CdgbkA8ekueXVka2aL4mgboGXcRn9geJw+nPC4tLIz\ngxtugEGDyh2Jc865Yssn+S3LUi2Z18OtJbWQ9AYwCxhvZu8A3c1sVlxkFqGbpLJ79llYvhz22qvc\nkTjnnCu2fBq8vCPpx0ArSZsDg4EX8tm4ma0EviupE/AvST/ImG+Scl7YGzZs2KrxVCpFKpXKZ7cN\ncuONod8+FaQ22Tnn3OqqqqqiqqqqKNvOpz+/DsBFwH5x0r+Ay+I1vfx3JP0G+Ab4OZAys5mS1iOU\nCHtnWb5kDV7mzIFNNoGPPoLOFXVHo3PONR0lbfBiZgvN7EJgb2AvM7son8QnaZ30zfGS2gH7EhrO\njAYGxsUGAo80NPhCGTkS9t/fE59zzlWKOqs9Je1I6L+vY3z9FXCSmU2oY9X1gOGSWhCS7F1mNk7S\nROABSScRb3VYjfhXm1m4t++KK8oZhXPOuVLKp9pzEnCamT0bX+8G3GRmfYsaWImqPV97DY49Ft57\nz6/3OedcY1bqLo2WpxMfgJk9BywvxM4bgxEjQvLzxOecc5UjZ8kv9t4OMABoB9wXXx8LLDazs4oa\nWAlKfmahocvo0dC3qOVY55xzq6uQJb/arvn9mZrP8xyaGG8Wz/Z8+WVo1w622abckTjnnCulnMnP\nzFIljKMs7r/fqzydc64S5dPgpQtwAtCT6mRZ9C6Nil3tuXIlbLQRPP009OlTtN0455wrkFJVe6Y9\nAbwIvEV4rFmzqPZ84QVYe21PfM45V4nySX5tzezsokdSYg88AMeU9Q5D55xz5ZJPtec5wHzgMWBJ\nerqZzSlqYEWs9ly2DDbcEJ57DjbfvCi7cM45V2ClrvZcDFxFeL5nujcHI89ujRqjJ56AXr088Tnn\nXKXKJ/n9GtjMzGYXO5hSuf12OPHEckfhnHOuXPKp9nwKOMLMFpYmpFX7LUq15/z5sMEGMGMGdOxY\n8M0755wrklJXey4C3pA0nuprfkW/1aFYxo6FXXf1xOecc5Usn+T3CN/udqjJ3uowZgwcemi5o3DO\nOVdOdVZ7rtbGpY2AO4H/IyTMW8zsOkldgRFAD2K3Rmb2Vca6Ba/2XLkS1lsPXnopPNPTOedc01HS\nXh0kTcsyfJDn9pcBZ5nZVsAuwOmS+gBDgLFm1gsYF18X3auvQrdunvicc67S5VPtuWNifA3gKGDt\nfDZuZjOBmXH8a0n/BTYA+gN7xsWGA1WUIAE++qhXeTrnnGtgtaek181s+3qu0xN4Btga+NjMusTp\nAuakXyeWL2i1pxn07g133w077lj38s455xqXkrb2jP36pbNQC6Af0LI+O5G0JvAQ8CszW6BENwpm\nZpKyZrlhw4atGk+lUqRSqfrstoZ33oHFi6FfvwZvwjnnXAlVVVVRVVVVlG3nc59fFdXJbzmhgcqf\nzGxqXjuQWgNjgH+a2TVx2hQgZWYzJa0HjDez3hnrFbTkd8klMG8e/OUvBdukc865EippyW91+vWL\nVZr/ACanE180GhgI/DH+zbyVouBGjoSbby72XpxzzjUFOUt+kgbG0eQCq7ozMrM769y4tBvwH0J3\nSOntXAC8AjwAbEwJbnV4911IpWD6dGhRZ/tW55xzjVGpSn478u2b2QUcCmxIuH+vVmb2HLlvp9gn\nnwAL4aGH4IgjPPE555wLciY/MxuUHpfUAvgRcD7wEvD74odWOCNHwp/+VO4onHPONRa1XvOLjVUG\nAucALwNH5dvQpbGYNg0++QR2373ckTjnnGssciY/SYOAwYQnsBxoZtNKFlUBjRoFhx0GrfK5nd85\n51xFqK3By0rgc+CLLLPNzPoWNbACNXjZdVf47W/hgAMKEJRzzrmyKWSDl9qSX8/aVjSzDwsRQC6F\nSH7Tp0PfvjBzJrRpU6DAnHPOlUVJWnsWO7mVwsMPh2d5euJzzjmX1Kwb/48aBUceWe4onHPONTZF\n7c9vdaxutecXX8Dmm8Nnn0G7dgUMzDnnXFmUtD+/purRR2H//T3xOeec+7YGJT9JlxQ6kEIbNQp+\n+MNyR+Gcc64xamjJb0JBoyiwuXPh+efhoIPKHYlzzrnGqEHJz8weK3QghfTww7DPPrDWWuWOxDnn\nXGOUT2e21xMecJ2+yGjAPGCCmT1axNgabMQIOOmkckfhnHOuscqnM9u/A1sADxIS4JHANKAr8IGZ\nnVmUwBrY2vOzz6B3b/j0U+jQoQiBOeecK4uSPOElsbOXge+b2fL4uhXwHLAbMMnM+tSy7m3AwcDn\nZrZNnNYVGAH0IEdffnG5BiW/c86BFSvg6qvrvapzzrlGrNS3OnQG1ky8XhPoGpPh4jrWvR3IfKrm\nEGCsmfUiPDR7SJ6x1mnJErjzThg0qO5lnXPOVa58+jq4Epgo6Zn4ek/gckkdgKdrW9HMns3yjND+\ncRsAw4EqCpQAH38cttoKNtusEFtzzjnXXOX1hBdJ6wM7ERq7TDCzGXnvICS/xxLVnnPNrEscFzAn\n/TpjvXpXex5/PKRScOqp9VrNOedcE1CSB1sndvYYcB/wqJktLMRO08zMJOXMcMOGDVs1nkqlSKVS\nObe1aBH8859w7bWFjNA551y5VFVVUVVVVZRt59PgJQUcCxwEvArcD4wxs7qu96XX70nNkt8UIGVm\nMyWtB4w3s95Z1qtXye+228JTXcaMyXsV55xzTUhJG7yYWZWZ/RLYDPgbcAyhk9uGGg0MjOMDgUdW\nY1sAmMENN3hDF+ecc/nJp8ELktoRGqocA2xPaKiSz3r3ERq3rCPpE+C3wB+ABySdRLzVof5h1/T8\n87BgAey33+puyTnnXCXIp9rzAWBn4ElCleczZray6IHVo9pz333h2GPh5z8vclDOOefKptQ3uR9A\nuC9vRXy9O3CcmZ1eiABq2W9eye+55+CEE2DqVGjdupgROeecK6eStvY0syclbS/peEIV5TTgoULs\nfHWZwZAhMHSoJz7nnHP5y5n8JG0BHE9o6fkF8dmeZpYqTWh1e/RRmD8ffvKTckfinHOuKclZ7Slp\nJTAGGGRmH8dp08xsk5IEVke159KlsPXWcN11cEDmA9Scc841O6W61eGHwDfAfyTdLGlvqrs1Krub\nbgqPMfPE55xzrr7yafCyJnAYoQr0B8CdwMNm9lRRA6ul5Pfll6HbomeegS23LGYUzjnnGouStvbM\n2HFX4ChCa8+9ChFALfvKmfxOPhnatg03tjvnnKsMZUt+pZQr+Y0ZE57k8tZb0LFjGQJzzjlXFiW9\n1aEx+eILOOUUuO8+T3zOOecarkmV/I47DjbcEP70pzIF5ZxzrmwqsuT3+OMwYULovcE555xbHU0i\n+T3/PPz0p/Dgg9C+fbmjcc4519TV2aVRuY0bB4cfDnfdBXvuWe5onHPONQdlS36SDpA0RdJ7ks7P\ntszbb4frfA89BPvvX+oIm6Zi9Xrc3Plxaxg/bg3jx638ypL8JLUEbgAOALYEjpfUJ3O5nXeGK6+E\nPfYodYRNl/9TNYwft4bx49YwftzKr1zX/HYC3jezDwEk3U94isx/kwtNnw5dupQ+OOecc81buao9\nNwA+SbyeHqfV4InPOedcMZTlPj9JRwIHmNnJ8fVPgJ3N7IzEMo3zBkTnnHNl09Tv85sBbJR4vRGh\n9LdKod6gc845l6lc1Z4TgM0l9ZTUhtBh7ugyxeKcc67ClKXkZ2bLJQ0C/gW0BP5hZv+tYzXnnHOu\nIBrtsz2dc865Yml0T3jJ5+b3SiVpI0njJb0j6W1Jg+P0rpLGSnpX0lOSOifWuSAeyymS9itf9OUn\nqaWkiZIei6/9uNVBUmdJIyX9V9JkSTv7catbPA7vSJok6V5Jbf24ZSfpNkmzJE1KTKv3sZK0Qzze\n70m6ts4dm1mjGQhVoO8DPYHWwBtAn3LH1VgGYF3gu3F8TWAq0Ae4EjgvTj8f+EMc3zIew9bxmL4P\ntCj3+yjj8TsbuAcYHV/7cav7mA0HfhbHWwGd/LjVecx6Ah8AbePrEcBAP245j9fuwHbApMS0+hyr\ndA3mK8BOcfwJwh0FOffb2Ep+q25+N7NlQPrmdweY2UwzeyOOf014KMAGQH/CSYr49/A4fhhwn5kt\ns/BAgfcJx7jiSNoQOAi4FUi3JPbjVgtJnYDdzew2CNfqzWweftzqMh9YBrSX1ApoD3yKH7eszOxZ\nYG7G5Pocq50lrQesZWavxOXuTKyTVWNLfnnd/O5AUk/Cr6WXge5mNivOmgV0j+PrU/MWkko+nlcD\n5wIrE9P8uNVuE+ALSbdLel3S3yV1wI9brcxsDvBn4GNC0vvKzMbix60+6nusMqfPoI5j2NiSn7e+\nyYOkNYGHgF+Z2YLkPAtl/tqOY8UdY0mHAJ+b2USqS301+HHLqhWwPXCTmW0PLASGJBfw4/ZtkjYD\nziRUy60PrBkf5LGKH7f85XGsGqSxJb86b36vdJJaExLfXWb2SJw8S9K6cf56wOdxeubx3DBOqzS7\nAv0lTQPuA/aSdBd+3OoyHZhuZq/G1yMJyXCmH7da9QNeMLMvzWw5MAr4Hn7c6qM+/5vT4/QNM6bX\negwbW/Lzm99rIUnAP4DJZnZNYtZowgV14t9HEtOPk9RG0ibA5oSLwhXFzC40s43MbBPgOODfZjYA\nP261MrOZwCeSesVJ+wDvAI/hx602U4BdJLWL/7P7AJPx41Yf9frfjN/V+bE1soABiXWyK3dLnywt\nfw4ktGJ8H7ig3PE0pgHYjXDN6g1gYhwOALoCTwPvAk8BnRPrXBiP5RRg/3K/h3IPwJ5Ut/b041b3\n8doWeBV4k1CC6eTHLa/jdh7hh8IkQoON1n7cch6r+wjXRpcS2nz8tCHHCtghHu/3gevq2q/f5O6c\nc67iNLZqT+ecc67oPPk555yrOJ78nHPOVRxPfs455yqOJz/nnHMVx5Ofc865iuPJz1UESStid0Zv\nS3pD0tnxZthi7e9xSR3rWGZgfHpFfbd9cezqZaqkf0vaMo91DpPUJ/G6StIO9d13Hvs5UdL1hd6u\nc4Xmyc9VikVmtp2ZbQ3sS3iYwtBi7czMDjaz+XUsdiLh2Y95kzQI2AXoa2ZbAFcAoyW1rWPVIwjd\nwawKMY99tapPbPlu17nGwJOfqzhm9gVwCjAIQg8Zkv4j6bU4fC9OHy5pVZdaku6R1F/SVpJejiXJ\nNyV9J3Mfkj6MHXL2VOgI9pZY6vyXpDUkHUV4BuQ9sceENSTtHcffkvSP+Ii/TOcBg8xscXwvY4EX\ngB/H/X6diOGo2CPD94BDgavi9jeNiwyI72GSpB3jOsMk3SXpOWC4pHUUOrN9JQ67xuV2kvRC3N7z\niUegJY/BwXGZrvX7hJwrPk9+riKZ2TSgpaRuhC5T9jWzHQjP/rwuLvYPQuks3bfd94DHgV8A15rZ\ndoRHKmV7+HqyBPQd4IZY6vwKONLMRhKeZfsjCz0mANwOHGNmfQk9KvwyucFYjdrBQj9mSROArbLs\n1+J7fZHwTMRzzGx7M/sgzm8X38NpwG2J9XoDe5vZj+OxuNrMdgKOIvSHCKEvyd1j7EOBy9NhxliP\nIHRCeqCFLn6ca1QaUq3hXHPTBrhB0rbACqAXgJn9R9JNktYhnPhHmtkKSS8AFyl0kDvKzN6vY/vT\nzOytOP4aoaubtPR1xy3icultDQdOB67NI36RX3Vj5jXO+yB0JiqpY0zwRnj26ZK4zD5An8Tl0bUk\ntQc6A3fGUq9R81yyF6FUu6+FTpeda3S85OcqUqz6WxGrQM8CPoslrn6EZJh2J+EJ8ScSS0dmdh+h\nGvEb4AlJP6hjd0sS4yuAlonXuZLWtxrjxGuIC+PT7JN2IDxEOXN77TI3UUec6fmLMuLYOV4v3c5C\n7xiLgMuAcWa2DeFYtEts43/AmoSE7lyj5MnPVZxY1XkzkG6V2BGYGcdPoGZyuoPQMamZ2ZS4/iZm\nNs3MrgceBbapbwjx74K4bwg9mfRU6AgVQsKtyrLuVcB1ktaIsewDfB+4N86fJam3pBaERi7phJbc\nVzqGY+M2diP0Nj6fbyfdp4DBq1YKpWPitj6N4z/N2O5HhJLynfm0RHWuHDz5uUrRLn2rAzAWeBK4\nNM67CRgo6Q1CaWVVVZ2ZfU7oi+32xLaOiY1XJhKutd2ZZX/fuvaW5fUdwM2SXo+vfwo8KOktYDkh\nQddcMSTcV4FJkqYAFwH9E9WUQ4AxwPNUJyeA+4FzY4OeTWMMi+O+bwJOSsSWjHcw0C827HkHODVO\nvxK4Iq7fMrGOhTBtKqERzoNZSqrOlZ13aeRcLeL1rbeA7cxsQbnjcc4Vhpf8nMshVilOJnSM6YnP\nuWbES37OOecqjpf8nHPOVRxPfs455yqOJz/nnHMVx5Ofc865iuPJzznnXMX5f52yaBZVtm5JAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cab5048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAADrCAYAAAAbp4NeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYlOW5x/HvTywoioJ6rNgVKyqW2KKrqMF4bInl2ILG\nWBK7JtZEUBNNLIlRj5oYCzaMGntBUVlAEwuKgiBqFASjgBxAEAUp9/njeRfGdcvs7szO7Ozvc11z\nMW+/39ll73nK+zyKCMzMzCrNEqUOwMzMrBic4MzMrCI5wZmZWUVygjMzs4rkBGdmZhXJCc7MzCqS\nE5yVFUlVkiaW8Pp3SJom6ZVWvu46kmZJUgHOdaGkWwsRVxOvu5qkoZJmSrq6Fa87XlKvVrrWnZIu\nb41rWcs5wbVzknaT9E9JMyT9n6SXJG1f6riaQ9Jxkoa14PjvA3sDa0bETvWcf0GWiGZJ+lDSKS2J\nuUZETIiIFaIAD6ZGxJURcWIh4mqik4ApEdE5In5V1w6SdpH0YpYEZ0h6XNJm+V6gngQT2as1tOa1\nrIWc4NoxSZ2BJ4E/A12AtYBLgbmljKuE1gXGR8ScBvZ5OUtEKwA/Bq6StE3rhFf21gXerW+jpJ2B\nZ4FHgDWA9YG3gZclrd8aAUpashCnKcA5rDVEhF/t9AVsD0xvYPtxwEvA1cA04COgd87244ExwEzg\nQ+CknG1VwCfAhcDnwDjgqJztPwRGZ8d+Apybc9xE4BxgMvApcFzOcSsCdwFTgPHAxaQ/OJsBXwPz\ngVnAtHruaU3gceD/gA+An2XrT6h1fN96Po9htda9CvxPzvJOwD+B6cBbwB4529YHhmb3PAj4X+Du\nbNt6wEJgiWy5Grgs+/xnkhLDyrX2/Qnwcfb5XpRznX51nLe+fZcF+mc/3zHAecDEBn4ndgFeB2YA\nrwE7Z+vvBL4hfTmaBexVx7HDgBvrWP800L+Bz3ghsCGphJh7jcey7eOAC0i/T9OA24Flav0engd8\nlt2rsv3/DUwF/g50ybneg9m+M4AhwOY52+4ALs/erwAMBq4r9f9lv+r5fS11AH6V8Ief/oNOzf44\n9c79T55tPy77g3JC9kfhFOA/Odt/CKyfvd8dmA1smy1XAfOAa4Clsu1fAhtn2z8Dds3er1jHcf2A\nDsB+2XlXzLbfRSoBdCKVGN4Dfppt61P7j2Md9zwUuBFYGtialCj3zOf42n98gR1JiWyjbHmt7PPs\nnS3vnS3XJKZ/AVcBSwK7Al8Ad2Xb1uO7Ce4DYCOgY/aH9Mpa+/4FWAboAcwBumfb+/LdBFffvr/P\nzr1iFv9IYEI99981u9+jSbU//0NKKF2y7XcAl9Vz7HKkLw971PO5flrXZ5ytWwhsUN81SF90Rmbx\ndyF9KahJQlWk36crSb+HHYEzSV9C1szW3QLcVyueTtm2PwEjcrbdQfrisTIpwdd5v36Vx8tVlO1Y\nRMwCdiO1KdwKTJH0mKT/ytnt44i4LdL/7ruANWq2R8TTETEuez8UeA74fq3L/CYi5mXbnwKOyNZ/\nA2whqXNEfBERI3KOmUf6w7EgIp4hJcbukjpkx18YEbMj4mPgWuDY7LgGq44kdSOVQM6PiG8i4m3g\nb6TSTaPHZ3aSNF3STOAVUoL6d7btGODpiBiYfSbPA8OB/SWtQyoxXxIR8yPiZVJJsr5rBnBHRPw7\nUpXpA0DtqtBLI2JuRIwkVfVt3cB91LfvYcAV2c/gP6Tq6vpi2h94LyLujYiFEXE/MBY4MGef+o7t\nSkqKn9WxbRKwSj3H1aX2NYJUMvxPREwHfgccmbN9IalEPi/7LE8Gfh0Rn0bEPFK1/KGSlgCIiDuz\n36+abVtLWiHnfGuRvoD8PSIuaULc1sqc4Nq5iBgbEcdHRDdgS9K32utydpmUs+9X2dvlASTtJ+mV\nrHPKdFKJbuWcY6dHxNc5yx9n54fUfvVDYLykakm5nTr+LyIW5ix/lV1zFdK36o9ztk0g/cHJx5qk\nqsvZzTwe4JWI6BIRnYHVgS0lXZFtWxc4LEuA07PPZNdsv5pr57bvNdZbdFLO+6/JPvd6tn9Vx/Z8\n9l2zVhyfNHCONUmfV67cn2lDppMSzRp1bFuDVHXaErn3MKFWTJ9HxDc5y+sBj+T8jMaQSperSeog\n6feS/i3pC1L1JyxOwCIl+o6kUrGVMSc4WyQi3iO1UWzZ2L6SlgH+Qapy+6+I6EJqS8n9dt1F0nI5\ny+sC/8muNTwiDgZWBR4llVAaM5VUulsvZ906LP6j3Fjvtk+BrpJyE0Hu8U0SEVOAh4EDslUTSFWD\nXXJeK0TEVaSSS1dJy9a6dql9BnTLWe5W346kn926tdYt+pk2JPtS8S/g8Do2Hw68kL2fTarOBEDS\n6rVPVc8l1qn1/tMGjplAqkbO/TktFxGfAUeRSqS9ImJFUrspLP69rqnteBZ4utbvt5UZJ7h2TFJ3\nSedIWitb7kaq2vlXHocvnb2mAgsl7QfsW8d+l0paStLupG++D2bLR0taMSIWkDoMLGjsgtm+DwC/\nk7S8pHWBs4F7sl0mA2tLWqqe4yeS2l6ulLSMpB7AT3OObxJJKwOHAO9kq+4BDpC0b1YS6Jg917dW\nVp06HOiX3f/OwH/TcFJujd56DwAXSlop+z04rYGYngY2kXSkpCUlHQFsSuqJC43HewHQR9LpklaQ\n1EXSb4HvkaoCIVWfbiFpa0kdSW2xuSYDG9RaJ+BUSWtJ6krqeHR/A3HcAlyRVRsjaVVJNdWsy5M6\nsUyT1Am4otaxAoiI00jtv09kcVoZcoJr32aR/ri8KulLUmIbCZybba/rmZ+ARe13Z5D+QE4jJcbH\nau07iVQ19Smp/e7kiHg/23YMMC6rBjqJ1HHhW9eox+mkb/kfkXrl3Utq+IdUChgNTJI0pZ7jjySV\nAD8llb4uiYgXG7jfXAHsXPMcHKlqa3IWExHxCXAQcBGp88oE0mdZ8//saGBnUg/Oy0m9976pdf7a\n18t9X3u5oTjz3fcyUgl2HKkN9cFaMS0+ScQ0UlI+l/TF5pfAf2fr67pu7eNfBn4A/Ij0+Y8ntQXu\nFhEfZvu8n8X0PCmBDKt1ztuAzbPqxYdzrntvFv+HpM45v23g/v9Mav98LmtL/RepwxCk39OPSaXS\nd7Jt9f0cTiJ9do9mNRpWZpT6DhThxOlbzRBSz62lSV16L8y+Yf2d7Jkj4PCImFGUIKxkJFWRqusa\nqvJq1yT9HRgTEZc2unMrkfRz0v/JPUsdi1lLFa0ElzWm7xkR25C6Ju8paTdSNcWgiNiE9I37gmLF\nYFZOJG0vaUNJS2RVugeS2h9LGdPqknbNYupOev7wkVLGZFYohXiqv145ve6WJj3TNJ30n3qPbH1/\nUndbJ7nK5CGNvm11UrXoyqRef6dkjyqU0tKkNqn1SQ82DwBuKmlEZgVStCpKgOy5kjdJoxDcHBHn\nSZqe9bhDkkhdp7sULQgzM2uXil2CWwhsI2lF4FlJe9baHpK+k2HrWmdmZhYRefcubpVelBHxBWkU\ni+2AyTXPtkhag9TbrK5jKvbVt2/fksfg+/P9tcf7q+R7aw/311RFS3CSVpG0UvZ+WWAfYASpe26f\nbLc+lLiR3czMKlMxqyjXAPpn7XBLkLqMvyBpBPCApBPIHhMoYgxmZtZOFS3BRcQooGcd66eRRllv\nt6qqqkodQlH5/tq2Sr6/Sr43qPz7a6qi9qJsLklRjnGZmVnpSCLKrZOJmZlZa3OCMzOziuQEZ2Zm\nFckJzszMKpITnJmZVSQnODMzq0hOcGZmVvZmNGPW0CYlOEkdJHVu+mXMzMyabs4cuOYa2Hjjph/b\naIKTNEBSZ0mdgFHAu5LOa/qlzMzM8rNgAfTvD927w0svwZAhTT9HPiW4zSNiJnAw8AywHnBs0y9l\nZmbWsAh4+mnYdlv461/h3nvh0Udh882bfq58xqJcUtJSpAT3vxExz/O1mZlZob32Gpx3HkyeDFde\nCQcdBMp7YK7vyqcE9xfSqP/LA0MlrQd80fxLmpmZLfbBB3DYYXDIIXD00TBqFBx8cMuSGzRjsGVJ\nAjpExPyWXbrBa3iwZTOzCjd5Mlx6KTzwAJxzDpx1Fiy3XP37N3Ww5UarKCV1AX5Canur2T+AM/K9\niJmZWY158+B//xd+9zs45hgYOxZWWaXw18mnDe5p4F/ASGAhIFKCMzMza5LBg+H002HNNWHYMNh0\n0+Jdq9EqSklvRsR3Ji4tJldRmplVlk8+gV/+El55Bf70p+a1sRVjPrj7JJ0kaQ1JXWteTQvLzMza\no9mz4bLLYJttYJNNYMyY1JmkpR1I8pFPFeUc4GrgYlIVJaQqyg2KFZSZmbVtCxbAHXdA376w++7w\n+uuw/vqtG0M+Ce5cYMOImFrsYMzMrG2LgGeeSc+zrbxyekh7hx1KE0s+Ce4D4OtiB2JmZm1XRBpO\n69JL4bPP4Kqr4IADWqcqsj75JLivgLckDQbmZusiIhp8TEBSN+Au4L9IVZp/jYjrJfUDfgZ8nu16\nYUQMbE7wZmZWWhEwaBBcfjlMmgQXXgjHHgtLLVXqyPLrRdmHbz8aIFKC69/IcasDq0fEW5KWB94g\nDfd1ODArIv7YwLHuRWlmVsYi4KmnUmKbNQsuvhiOOAKWzKfY1EwFfdBb0pLA8RFR1dRAImISMCl7\n/6Wkd4G1ak7d1POZmVl5GDsWTjkFpk2D3/wGfvxjWKIMZxdtMKRsOK4FklZqyUWy8Su3BV7JVp0u\n6W1Jt7X03GZm1jrmzIF+/WC33eBHP4IRI9IYkuWY3CC/NrjZwChJg7L3kEcbXI2sevIh4MysJHcz\ncFm2+XLgWuCE2sf169dv0fuqqiqqqqryuZyZmRXB4MGp1Lb55imxdetW/GtWV1dTXV3d7OPzaYM7\nro7VjbbBZccuBTwJPBMR19WxfT3giYjYqtZ6t8GZmZWBSZPgggvgxRfh+uvTCCSlUvDBliPizmYG\nIuA2YExucpO0RkR8li0eQpol3MzMysicOXDddXDNNXD88TB6NKywQqmjapp8ZhPYBLgC2BxYNlsd\nEdHYSCa7AscAIyWNyNZdBBwpaRtSr8xxwMnNCdzMzAovAh5+GH71K+jRI40dudFGpY6qefJpg7sD\n6Av8EagCjgc6NHZQRLxE3Z1YnmlCfGZm1kpqZtSeNg1uvRV69Sp1RC2TT9+XZSPieVJ73ccR0Q/Y\nv7hhmZlZa1i4EJ58EvbYI/WIPPJIePPNtp/cIM/BliV1AP4t6TTgU6BTccMyM7NimjsX7rkHrr0W\nOnZMVZKHHVbcB7VbWz69KHcE3gVWInXr7wxcFRGvNHhgS4JyL0ozs6KYPh1uuQVuuAG23jrN0bbX\nXqUdMzJfTe1FWW+Ck/RcROybvb8oIq4oUIyNB+UEZ2ZWUB9/nCYaveuuNAjyueemTiRtSSEnPF01\n5/1hzQ/JzMxKZdSo1K7Ws2caAHnkSOjfv+0lt+Yo0wFWzMysJebPh9/+NlU/9uwJH30EV18Na69d\n6shaT0PNiRtIepw0MPL6kp7I2RYRcWBxQzMzs+b44AP4yU+gU6fUI7I1htUqRw21wVU1cFxExJCi\nRITb4MzMmiMCbr4ZLrkE+vaFU08t34GQm6NgQ3VFRHVBIjIzs6J7/30444z0kPZLL8Gmm5Y6otKr\noNxuZtb+vPNO6kSy665QVQUvv+zkVsMJzsysDRoxIk002qtXep7tww/TqP9LLVXqyMqHE5yZWRvy\n6adpstH990+lto8+Somtc+dSR1Z+6m2Dq91rktSbctGye1GambWeCBgwAM4+O008eu+9sOyyjR/X\nnjX0mMC12b+HAKsD95CS3JHA5CLHZWZmmSlT4Oc/h7Fj4emnYbvtSh1R25DPWJRvRMR2ja0raFB+\nTMDMDEhzs516anqu7dJL08DI7VXBZ/QGlpO0YUR8mF1gA2C55gZoZmaN+/xzOPNMeP11eOih1N5m\nTZNPJ5OzgcGShkgaAgwGzipuWGZm7VNEal/bcktYay14+20nt+ZqtIoSQFJHoHu2ODYi5hY1KFdR\nmlk7NHFi6kAycSLcdhvssEOpIyovhZxNoOaEnYBfAadFxNvAOpL+uwUxmplZjgUL4Kab0qDIO+0E\nw4c7uRVCPm1wdwBvALtky58CDwFPFisoM7P2oroazjorPcc2ZAhsvnmpI6oc+bTBbRgRfwC+AYiI\n2fmcWFI3SYMljZb0jqQzsvVdJQ2S9L6k5ySt1Pzwzczapo8+Sg9sH3ccXHyxk1sx5JPg5kpa9Dih\npA2BfNrg5gFnR8QWwE7AqZI2Ay4ABkXEJsAL2bKZWbswcyacf36qgtx+e3j3XTjsMFDeLUuWr3wS\nXD9gILC2pPuAF4HzGzsoIiZFxFvZ+y+Bd4G1gAOB/tlu/YGDmx62mVnbMmcOXHcdbLwxTJ6cZtq+\n6CKPRlJM+faiXIVUCgN4NSI+b9JFpPWAIcCWwISI6JKtFzCtZjlnf/eiNLOKMH8+3H039OuXBkX+\n3e9gq61KHVXbVPAHvSW9EBG9yOlUkrMun4CWB/4BnBkRs5RTDo+IkFRnJuvXr9+i91VVVVRVVeVz\nOTOzshABjz6a2tdWWQXuu8/PszVVdXU11dXVzT6+oRm9lyWNWDIYqMrZ1BkYGBGNzjgkaSlSYnwm\nIq7L1o0FqiJikqQ1gMG1z+USnJm1ZePGpbEj//MfuOoq6N3bbWyFUMjn4E4GhpMe8H4j5/U4cGMe\ngQi4DRhTk9wyjwN9svd9gEfzDdbMrJzNnw/XXps6kOy5J7z5Juy3n5NbqeQz2PLpEXFDk08s7QYM\nBUaSptsBuBB4DXgAWAcYDxweETNqHesSnJm1KSNGwIknwoorwl/+AhttVOqIKk8xBlsOSV0iYnp2\ngS7AkRFxU4MHRbxE/SXEvfMN0MysnE2bljqO3HMP/OEP0KePS2zlIp/HBE6sSW4A2fuTiheSmVn5\n+/rr1L7WvTvMnp26/R93nJNbOcmnBLeEpCUiYiGApA7AUsUNy8ysPC1YAHfdBZdcktrahg2DTRvt\ncmelkE+Cexa4X9JfSDN6n0x68NvMrF0ZNiz1juzSBR54AHbeudQRWUPy6WTSgVQlWfPc2yDgbxGx\noGhBuZOJmZWROXPgN79J87TdeCMccoirIkuh4J1MImKBpP6k59XGtig6M7M25q234NhjYZNN0uSj\nq65a6ogsX/nMB3cgMIKsWlLStpIeL3ZgZmalNH8+XHEF7LMPnHcePPSQk1tbk08bXD/ge6QRTYiI\nEZI2KGZQZmal9NprcMYZ0KkTvPEGrLNOqSOy5sjnMYF5tR/EBhYWIxgzs1KaMAGOOSa1sZ14Igwa\n5OTWluWT4EZLOhpYUtLGkm4A/lnkuMzMWs2sWfDrX8O228KGG8J778EJJ8AS+fyFtLKVz4/vNGAL\n0iSnA4CZwFnFDMrMrDUsXAi3354e1p44MXUiufRSWH75UkdmhdDQbAJ3R8Sxks6qNVhy8YPyYwJm\nVmTDh8Opp6ZS2g03pNm1rbwVcjaB7SStCfxUUtfar5aHambW+qZOhZNOggMOSA9tv/yyk1ulaqgX\n5S3AC8AGpGlyckW23sysTZg3D/7611QFedRR8O67sNJKpY7KiimfkUxuiYhTWimemmu6itLMCmL+\n/DQCyWWXwQYbwB//CFttVeqorDmaWkXZaILLTtoBWI2cEl9ETGhWhPkE5QRnZi20YAHcf38qsa2x\nRkpwe+xR6qisJQo+VJek04G+wBQgd/xJfwcys7ITAf/4Rxrtv0sXuOWWNLu2x45sf/KpovwQ2DEi\n/q91QnIJzsyap7o6Dau1YEEaZmvffZ3YKkkxZvSeQHr2zcysLI0cCRdcAGPHptm1jzjCD2lbfglu\nHDBY0lPAN9m6iIg/Fi8sM7PGTZyYprF55hm4+GJ49FFYeulSR2XlIt8S3ARg6ewl0mMCZmYl8eWX\n8Ic/wE03pWfZPvgAOncudVRWbvKZD65fK8RhZtaoBQvgzjtTqa1XrzRXW7dupY7KylW9CU7SEw0c\nFxFxYGMnl3Q7sD8wJSK2ytb1A34GfJ7tdmFEDMw7YjNrl154Ac45B1ZcER57DHbYodQRWblraCzK\nqgaOi4gY0ujJpe8DXwJ35SS4vsCshtrw3IvSzGqMGwfnnptKa1dfDT/6kXtGtlcF60UZEdUtDSYi\nhklar45N/vU0swbNng2//z3cfHMqud13H3TsWOqorC0pVUfa0yW9Lek2SR4NzswWiYABA2DTTeGj\nj1LJ7aKLnNys6fLpRVloNwOXZe8vB64FTqi9U79+/Ra9r6qqoqqqqhVCM7NSGjkSTjst9ZIcMAB2\n263UEVkpVVdXU11d3ezjiz4fXFZF+URNG1w+29wGZ9a+zJiRhtb6+9/TmJE/+xl06FDqqKzclP18\ncJLWyFk8BBjV3HOZWdtWM6P2ppum6WzGjIGTT3Zys8Io6nxwkgYAewCrSJpIGrS5StI22TnGASc3\nI24za+PeeCPNqC3B009Dz56ljsgqjeeDM7NWNX06/PrXacT/K6+EPn08bqTlp5BVlABExCmStpZ0\nuqTTJG3dshDNrD1auDCNQrL55qmn5JgxcPzxTm5WPPnMB3cmcCLwMOn5tXsk3RoR1xc7ODOrDG+9\nlXpHfvMNPPEEbL99qSOy9iCfKspRwE4RMTtb7gS8UlevyIIF5SpKs4owfnyqjnz++TSztntHWksU\nvIoys7Ce92Zm3zF1Kpx9Nmy3HWy0URrt370jrbXl86D3HcCrkmqqKA8Gbi9qVGbWJs2eDddfD9de\nmyYdHT0aVl+91FFZe9VoFSWApO2A3Uhd+4dFxIiiBuUqSrM2Ze5c+MtfUq/I3XeH3/4WNt641FFZ\npSnYYMu5IuINvvssnJm1c/Pnp56Rl18OPXqkmbW32abUUZklpRiL0szauIUL07Bal1wC66wD998P\nO+9c6qjMvs0JzszyFpFKaTWj+99yS5pZ26wcNZjgJC0JDIqIPVspHjMrUy+/DBdemHpIXnEFHHSQ\nJx618tbgYwIRMR9Y6DnbzNqvsWPhwAPhqKPgpz+FUaPg4IOd3Kz85VNFORsYJWlQ9h4gIuKM4oVl\nZqU2a1aauubOO1PJ7cEHYZllSh2VWf7ySXAPZ6+afvvKeW9mFSYC7rsPzjsP9t0X3nkHVlut1FGZ\nNV2+z8EtB6wTEWOLH5KfgzMrlREj4Mwz0wPbN97onpFWXgo+VJekA4ERwMBseVtJjzc/RDMrJxMn\nwtVXp/nY9t8/tbW99pqTm7V9+Qy2/CawFzA4IrbN1r0TEVsWLSiX4MyKasaM9OzagAGpCvKQQ1Ji\n22MPjxdp5asYI5nMi4gZ+naXKQ+4bNYGTZgA112XOo7svTeccw707u3OI1aZ8klwoyUdDSwpaWPg\nDOCfxQ3LzAppxAi45hoYODBNMvr229CtW6mjMiuufKbLOR3YApgLDABmAmcVMygzK4whQ2CffeCA\nA9IYkR99lBKdk5u1B3n1ogSQtCLp+beZxQ3JbXBmLTVkCPTrl6okL74YjjkGll661FGZtUzB2+Ak\n7UCa/61ztjwDOCEihjc7SjMriurqlNg++STNpH300bDUUqWOyqw08qmivB34RUSsGxHrAqeS54Sn\nkm6XNFnSqJx1XSUNkvS+pOc8DJhZyw0dClVV8LOfpTa2sWPhuOOc3Kx9yyfBzY+IYTULEfESMD/P\n898B9K617gLSAM6bAC9ky2bWDP/6V2pjO+649Bo7Fvr0gSU9T4hZ/W1w2SzeAMcCy5I6mAAcAcyJ\niLPzuoC0HvBERGyVLY8F9oiIyZJWB6ojYtNax7gNzqwBr70GffvCmDGpKtKlNWsPCtkGdy3fHn+y\nb877lmSf1SJicvZ+MuBR7szyNHUq/PKXMGhQ6jzy6KN+hs2sPvUmuIioKvbFIyIk1Zks+/Xrt+h9\nVVUVVVVFD8esbEVA//5w/vmp48jYsbDCCqWOyqy4qqurqa6ubvbx+QzV1QX4CbAeixNi3tPl1FNF\nWRURkyStQRoCzFWUZvV47z045RSYORP++lfYbrvGjzGrRAUfbBl4GlgXGAkMB97IXs31ONAne98H\neLQF5zKrWDNnwm9+A7vummbPfvVVJzezpsinr9UyEXFOc04uaQCwB7CKpInAJcDvgQcknQCMBw5v\nzrnNKtVXX6Wpaq65Bn7wA3jzTVhnnVJHZdb25FNF+UvS8FxPkIbrAiAiphUtKFdRWjs0dy7ceitc\ncQXsskuaTXvzzUsdlVn5KMZsAnOAq4GLWTyLQAAbND08M6tt6lS4+274859hs83giSdcFWlWCPmU\n4MYBO0TE1NYJySU4q3wLF8LgwanENnBgGgz55z9PJTczq1sxSnAfAF83PyQzqzFlCtx2G/ztb6mb\n/4knws03Q5cupY7MrPLkk+C+At6SNJjFbXB5PyZgZvD663DDDan68cc/TrNpb789KO/vombWVPkk\nuEf5bld+1x+aNWLuXHjwwdQjcvJk+MUv0mzaXbuWOjKz9iHv+eBak9vgrC176y244w647740yejp\np8P++0OHDqWOzKxtK8Z8cOPqWB0R4V6UZplp01JCu/321Cvy+OPTgMjrr1/qyMzar3yqKHfIed8R\nOBRYuTjhmLUdCxbA88+npDZwIPzwh/CHP8Bee7m0ZlYOmlVFKenNiOhZhHhqzu8qSitbH34Id96Z\nXqutBj/9KRx5pHtCmhVbMaoot2Nxp5IlgO0Bfz+1dmXqVHjooVQN+e67aUT/J5+ErbcudWRmVp98\nHvSuZnGCm08aP/KaiHivaEG5BGdl4Msv4fHHU1IbNgz22y+V1Hr39hxsZqXQ1BKce1Ga5fjmm9Se\nNmAAPPNMGsn/qKPgwAM9/5pZqRUswUmqmdImd4dFs3lHxF3NDbLRoJzgrBUtWABDh6aS2iOPwBZb\npJLaoYfCKquUOjozq1HINrgd+O4D3QIOANYGipbgzFrD22/DPfekxLbaaimpjRgB3bqVOjIzK4S8\nqiglLQEcBZwPjAF+FxEjixaUS3BWJJ98khLaPffAF1/AMcekDiOelsas/BW0F6WkpUizbv8SeBU4\ntJidS8zRqK/LAAAPC0lEQVSKYfz4VPX48MMwenSqerzxRthtN1ginzntzaxNaqgN7jTgDOAF4KqI\nqGtEk+IE5RKctUBESmSPPJJen3ySOokccgj06gUdO5Y6QjNrjkJ2MlkITAE+r2NzRESP5oWYR1BO\ncNZE8+fDyy/DY4+lrv3z5qWEdsghqSfkkvmM2WNmZa2QVZQea9LK2tdfpy79jzwCTz8N664LBx2U\nHsjeemtPRWPW3vk5OGtTZs9Oyeyhh1Jy2357+NGPUhWkez+aVbY286C3pPHATGABMC8idszZ5gRn\ni4wfnwY1fuaZ9O9OO6WOIgcfDKuuWurozKy1tKUENw7YLiKm1bHNCa4dmz49JbIXXkj/zpoFe+8N\n++wDBxwAK3suC7N2qeCDLReZW0kMgIULYfBg+NvfFg+RtffecOqpsOWWbk8zs6Zr7nQ5l0ZE3xZd\nWPoI+IJURfmXiLg1Z5tLcO3Ep5+maWduuw2WXx5OPDE9eO2pZ8ysttYqwQ1v5nG5do2IzyStCgyS\nNDYihtVs7Nev36Idq6qqqKqqKsAlrRzMnAmPPpoGNH71VTjsMLj//tRhxCU1M6tRXV1NdXV1s48v\ni16UkvoCX0bEtdmyS3AVZs6cVPU4YAA8+yzssUcapf+AA6BTp1JHZ2ZtQTEmPL2BNOhyzUmDVLU4\nPCIea2aQywEdImKWpE7AvsClzTmXla+FC+Gll+Duu+Ef/0jPph11FNxyC3TtWurozKzS5VNF2RHo\nDjxISnI/BsYBW0vaMyLOasZ1VwMeUaqPWhK4NyKea8Z5rAy9+25KavfeC507w7HHwsiRsPbapY7M\nzNqTfGb0fpXUXjY/W14SeAnYDRgVEZsVPChXUbYZX3+dpph59VV47bX079y5qaR27LHQo2gDuplZ\ne1OMTiYrAcsDM7Ll5YGuETFf0pxmxGht2Oefw7BhaYLQYcNg7Ng01cyOO0Lv3nDJJdC9u0fpN7PS\nyyfBXQWMkDQkW94DuCJrO3u+aJFZWZg6FV58MT2jNnRoGpl/111TJ5EbboCePT06v5mVp3wnPF0T\n2JHUwWR4RPynqEG5irJkvvoqlcyefz69PvooJbOqqvTvNttAhw6ljtLM2qOCD9Ul6QlgAPBYRMxu\nYXz5BeUE12rmz4fhwxcPizV8eCqV9eqVRhLZYQdYaqlSR2lmVpwEVwUcAfwQeB24H3gyIorW/uYE\nV1wTJsCTT8Jzz0F1dZpmZu+90+v7308jipiZlZuiDbac9Z7cEzgR6B0RnZsXYl7XcoIroAh46600\nEehjj8HHH8P++6dOIb16wWqrlTpCM7PGFWWoLknLAgcChwM9gf7NC89ay8cfp04hQ4bAoEFpRuuD\nDoI//ckzXJtZ+5BPFeUDwPeAgaTqySERsbCoQbkE12RTp8ITT6QqxyFD0vNpu++eXr16wWabeZxH\nM2vbitEG1xsYFBELsuXvA/8TEae2KNKGr+kEl4cpU+CRR+DBB+H119N8afvsk3o7du/uhGZmlaXg\nVZQRMVBST0lHkqooxwH/aEGM1kzffANvvpnGd3zqqTSCyH77wS9+kdrTlluu1BGamZWPehOcpO7A\nkaQelJ+TjUUZEVWtE1r7FpHmSnv7bXj55fQaPhw23ji1oZ15JvzgB7DssqWO1MysPNVbRSlpIfAk\ncFpETMjWjYuI9YseVDuqovzmm9Rtf9y4NOzVO+/A6NHptfTSsNVWsMsuKanttBOsuGKpIzYzK41C\nVlH+iFSCGyppIItnE7Bm+OqrVL346qupVDZuHIwfn9rR1loL1lsvtZv16AFHHglbbAGrrlrqqM3M\n2q58OpksDxxESnZ7AncBjxRzepu2XIKLgMmT4f334b33UrXia6+l5S22SIMS9+wJG2yQktraa7vL\nvplZPor2oHd28q7AoaRelHs1I758r1P2CS4iPWs2cmR6jR4NH3yQEtkyy8Amm6RXz54pqW2zTVpv\nZmbNU9QE11rKMcFNmZKeMRs6NI0KMmoUrLBCqlLs0SNNGdO9e+oE4tmqzcwKzwmuQKZPTwmtZqqY\nTz5J4zRWVaVSWY8esPLKJQ3RzKxdcYJrpvnzUweQ556DZ59NVY677gp77gl77QXbbuu2MjOzUqrY\nBDdwIDzwQHGuN21aGt5q3XXTs2X77gu77eY2MzOzclKxCW7UqDQcVTF06pSGt1p99eKc38zMWq5i\nE5yZmbVvTU1wSxQzmPpI6i1prKQPJJ1fihhKqbq6utQhFJXvr22r5Pur5HuDyr+/pmr1BCepA3Aj\n0BvYHDhS0matHUcpVfovoe+vbavk+6vke4PKv7+mKkUJbkfg3xExPiLmkeaYO6gEcZiZWQUrRYJb\nC5iYs/xJts7MzKxgWr2TiaQfA70j4sRs+RjgexFxes4+7mFiZmbfUdAJT4vgP0C3nOVupFLcIk25\nATMzs7qUoopyOLCxpPUkLU2aUPXxEsRhZmYVrNVLcBExX9JpwLNAB+C2iHi3teMwM7PKVpYPepuZ\nmbVUSR70bkilPQQu6XZJkyWNylnXVdIgSe9Lek7SSqWMsbkkdZM0WNJoSe9IOiNbXyn311HSq5Le\nkjRG0pXZ+oq4vxqSOkgaIemJbLli7k/SeEkjs/t7LVtXSfe3kqSHJL2b/Y5+rxLuT1L37GdW8/pC\n0hlNvbeySnAV+hD4HaT7yXUBMCgiNgFeyJbbonnA2RGxBbATcGr286qI+4uIOcCeEbEN0APYU9Ju\nVMj95TgTGAPUVOdU0v0FUBUR20bEjtm6Srq/PwNPR8RmpN/RsVTA/UXEe9nPbFtgO+Ar4BGaem8R\nUTYvYGdgYM7yBcAFpY6rAPe1HjAqZ3kssFr2fnVgbKljLNB9PgrsXYn3BywHvA5sUUn3B6wNPA/s\nCTyRrauk+xsHrFxrXUXcH7Ai8FEd6yvi/nLuZ19gWHPuraxKcLSfh8BXi4jJ2fvJwGqlDKYQJK0H\nbAu8SgXdn6QlJL1Fuo/BETGaCro/4E/Ar4CFOesq6f4CeF7ScEknZusq5f7WBz6XdIekNyXdKqkT\nlXN/Nf4HGJC9b9K9lVuCa3c9XiJ9FWnT9y1peeAfwJkRMSt3W1u/v4hYGKmKcm1gd0l71treZu9P\n0n8DUyJiBFDns6dt+f4yu0aq5tqPVIX+/dyNbfz+lgR6AjdFRE9gNrWq7Nr4/ZE9SnYA8GDtbfnc\nW7kluEYfAq8QkyWtDiBpDWBKieNpNklLkZLb3RHxaLa6Yu6vRkR8ATxFag+olPvbBThQ0jjSN+S9\nJN1N5dwfEfFZ9u/npDacHamc+/sE+CQiambKfIiU8CZVyP1B+mLyRvbzgyb+7MotwbWXh8AfB/pk\n7/uQ2q7aHEkCbgPGRMR1OZsq5f5WqemlJWlZYB9gBBVyfxFxUUR0i4j1SdVAL0bEsVTI/UlaTtIK\n2ftOpLacUVTI/UXEJGCipE2yVXsDo4EnqID7yxzJ4upJaOLPruyeg5O0H3Adix8Cv7LEIbWIpAHA\nHsAqpDrjS4DHgAeAdYDxwOERMaNUMTZX1qNwKDCSxVUFFwKvURn3txXQn/RFcAlSKfVqSV2pgPvL\nJWkP4NyIOLBS7k/S+qRSG6TqvHsj4spKuT8ASVsDfwOWBj4Ejif97Wzz95d9KfkYWL+m6aOpP7uy\nS3BmZmaFUG5VlGZmZgXhBGdmZhXJCc7MzCqSE5yZmVUkJzgzM6tITnBmZlaRnOCsXZC0IJt2451s\n+ptzsgfVi3W9pyR1bmSfPtloDE0996+z6ULek/SipM3zOOag3Jk5JFVL2q6p187jOsdJuqHQ5zVr\nDic4ay++ijT9xpakEUn2A/oW62IRsX9EzGxkt+OANZtyXkmnkaYm6hER3YErgcclLdPIoYeQpqBa\nFGIe11qyKbHle16z1uIEZ+1ONq7dScBpkGZCkDRU0hvZa+dsfX9JB9UcJ+leSQdK2kJpItQRkt6W\ntFHtayhNtNk1O/e7kv6alR6fVZpI9VBge+DebCT4jpJ6Ze9HSrotG66utvOA0yLNVUdEDAL+CRyd\nXffLnBgOzUaa35k0YO3V2fk3yHY5NruHUZJ2yI7pJ+luSS8B/bPhyh6S9Fr22iXbb0dJ/8zO93LO\ncFG5n8H+2T5dm/YTMisMJzhrlyJiHNBB0qqkIdT2iYjtSGMyXp/tdhuplIWkFUnzFT4FnAL8ORZP\nxljXgOC5JZmNgBuz0uMM4McR8RBp7NWjspHgIU2Oe3hE9CANLfXz3BNmVZ6dImJ8rWsNJ81TV/u6\nkd3rv0hj+P0yInpGxEfZ9mWze/gFcHvOcZsCvSLi6Oyz+FOkyUIPJQ0LBfAu8P0s9r7AFTVhZrEe\nApwP7BcR0+r4fMyKrjlVEGaVZmngxmxcvwXAJgARMVTSTZJWIf1xfygiFkj6J3CxpLWBhyPi342c\nf1xEjMzev0GaALdGTTtg92y/mnP1B04lzdjcGJFf1WDtNscBABExTFLnLIkH8HhEzM322RvYLKe5\ncgVJywErAXdlpdfg239L9iKVTveJiC8xKxGX4KxdyqrpFmTVlWcDn2Ulp+1JCa/GXcCxpJLc7QAR\nMYBU5fc18LRqzRFXh7k57xeQBsOtUV9i+k4HmKxNb3Y2iHCu7UijyNc+37K1T9FInDXbv6oVx/ey\n9stts9kHvgIuB16IiK1In8WyOef4EFielLTNSsYJztqdrFryFqCmt19nYFL2/id8OwHdCZxFml9x\nbHb8+hExLiJuIM0MsVVTQ8j+nZVdG+A9YD1JG2bLxwLVdRx7NXC9pI5ZLHsDuwL3ZdsnS9pU0hKk\njiU1SSv3WjUxHJGdYzdgRpZAayfW54AzFh2USrlk5/o0e398rfN+TCrx3pVPD0+zYnGCs/Zi2ZrH\nBIBBwEDgsmzbTUAfSW+RSh2LqtUiYgowhtQ+VuPwrMPICFLb1111XO87bWF1LN8J3CLpzWz5eOBB\nSSOB+aQk/O0DU1J9HRglaSxwMXBgTpXiBcCTwMssTkAA9wO/yjrRbJDFMCe79k3ACTmx5cZ7BrB9\n1plmNHBytv4q4Mrs+A45x0QKM94jdXx5sI4Sp1mr8HQ5Zg3I2ptGAtvWzEllZm2DS3Bm9ciq/8YA\n1zu5mbU9LsGZmVlFcgnOzMwqkhOcmZlVJCc4MzOrSE5wZmZWkZzgzMysIv0/lS2uxmR+fIkAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1090a45f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(1,figsize =(7,7) )\n",
    "plt.subplot(211)\n",
    "plt.plot(list(range(1,1001)),num_inf/100)\n",
    "plt.xlabel('Days into Outbreak')\n",
    "plt.ylabel('Avg. Number of Infected Farms')\n",
    "plt.title(\"Outbreak: Ideally All Farms Infected by 1000 Time Steps\")\n",
    "\n",
    "plt.figure(2,figsize = (7,7))\n",
    "plt.subplot(211)\n",
    "plt.plot(list(range(1,71)),num_inf[0:70]/100)\n",
    "plt.xlabel('Days into Outbreak')\n",
    "plt.ylabel('Avg. Number of Infected Farms')\n",
    "plt.title(\"Snapshot of Beginning of Outbreak\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}