{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Python using scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Author*: Simon Frost\n", "\n", "*Date*: 2018-07-12" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from scipy.integrate import ode, solve_ivp" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def sir_ode(times,init,parms):\n", " b, g = parms\n", " S,I,R = init\n", " # ODEs\n", " dS = -b*S*I\n", " dI = b*S*I-g*I\n", " dR = g*I\n", " return [dS,dI,dR]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "parms = [0.1,0.05]\n", "init = [0.99,0.01,0]\n", "times = np.linspace(0,200,2001)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "sir_sol = solve_ivp(fun=lambda t, y: sir_ode(t, y, parms), t_span=[min(times),max(times)], y0=init, t_eval=times)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sir_out = pd.DataFrame({\"t\":sir_sol[\"t\"],\"S\":sir_sol[\"y\"][0],\"I\":sir_sol[\"y\"][1],\"R\":sir_sol[\"y\"][2]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Visualisation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "plt.style.use(\"ggplot\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAENCAYAAABUyo/yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8U/X+x/HXSdJ0r3RQRqFQQDalspcIBRfI+ImgKCB4FZCtF714FRRR5KKgCAKKgBNQGQIqCA72LGXIEAq10AFt092mzTi/PyJVBKSFtknbz/PxyKNNcpLzTsW+e9b3q6iqqiKEEEKIW6ZxdAAhhBCiopMyFUIIIW6TlKkQQghxm6RMhRBCiNskZSqEEELcJilTIYQQ4jZJmQohhBC3ScpUCCGEuE1SpkIIIcRtkjIVQgghbpPO0QFuV2Ji4i29LjAwkNTU1FJOc/skV8lIrpJz1mySq2RuJ1eNGjVKOY2QLVMhhBDiNkmZCiGEELdJylQIIYS4TeVyzHThwoVER0fj6+vLW2+9dc3zqqqybNkyDh8+jKurK2PGjKFevXrlEU0IIYS4beWyZdqtWzemTp16w+cPHz5McnIy7777Lk899RQffvhhecQSQgghSkW5lGmTJk3w8vK64fMHDx6ka9euKIpCw4YNyc3NJT09vTyiCSGEELfNKS6NMRqNBAYGFt0PCAjAaDTi7+9/zbJbt25l69atAMyaNeuq15WETqe75deWJclVMpKr5Jw1m+QqGWfNVVU5RZmWRFRUFFFRUUX3b+U6K+/XX0f366+Ya9XCUrculnr1sNati6V2bXB1Lc24JVYZr2krS5Kr5Jw1m+QqGbnO1Lk4RZkaDIar/lGkpaVhMBjKbH36gwfR7NuH598eVzUarLVqYalXD0ujRpgbN8bcuDGW+vUdXrJCCCGcl1OUaevWrfn+++/p1KkTZ86cwcPD47q7eEtLxrx5GC5fJu/IEbTnz6M7dw7d+fNoL15EFx+PLj4efv65aHlVp8NSv769XJs1wxwZibl5c1R39zLLKIQQouIolzKdN28eJ06cIDs7m1GjRvHwww9jsVgA6NWrF61atSI6Oprx48ej1+sZM2ZMmeax1q6NGhlJbuvWVz9RUGAv09hYdCdP4nLiBC4nT6KNi8Pl1ClcTp2CtWsBe8GaGzfGHBlJYWQkhW3aYK1Tp0xzCyGEcE6Kqqqqo0PcjvIYm1fJy0N3+rS9XI8cQR8dje70aRSb7arlLKGhFHTqRGGnThR07IgtJKRMc5UnyVUyzpoLnDeb5CoZOWbqXJxiN6+zUz08MLdqhblVKxgyBAAlN7eoWF2io3Hdtw/dhQvoVq7Ec+VKAMz161PQvTumnj0pbNsWdPLjFkKIykh+u98i1dOTwo4dKezY0f6A1Yru5Elcd+7Eddcu9Hv34nL2LC5nz+K1ZAk2Pz9Md9+NqWdPCrp3R/X2duwHEEIIUWqkTEuLVoulWTMszZqRO2oUmM3oo6Nx++EHXH/4AZezZ/FYuxaPtWtR3dww9ehBfr9+mLp3Bzc3R6cXQghxG6RMy4qLC4Xt2lHYrh38979oz53D7YcfcNu8Gdd9+3DftAn3TZuweXtjuvde8gYNorB9e0enFkIIcQtk1phyYq1Xj9ynnyZtzRqS9+8n86WXKGzeHE12Nh5ffkngQw8R3LUrmrlz0RiNjo4rhBCiBKRMHcBWsya5o0aR+v33XPrlF7InTMAaEmK/3vWFF6h25534jRmDS3S0o6MKIYQoBilTB7PWr0/2lClc2rePtGXLsN13H1gseKxfT1CfPgT064fbd9+B1eroqEIIIW5AytRZ6HQU9OqFZd06Lu3dS/bYsdh8fXE9cADDk08S3LUrHp98AoWFjk4qhBDib6RMnZCtZk2y//MfLh04QOarr2IJDUUXF4ffCy8Q3KULHp9/Dmazo2MKIYT4g5SpE1M9PckdOZLLO3diXLgQc8OG6C5exO/f/ya4a1fcV62S3b9CCOEEpEwrAp0OU9++pGzdSvqCBZjDw9HFx+M/eTJB99yD6/btjk4ohBBVmpRpRaLVkt+vHyk//UT6u+9iCQ3F5eRJAh55BMPjj6M7c8bRCYUQokqSMq2ItFry/+//uPzzz2RNnYrNywu3H38kqEcPfF5+GSU729EJhRCiSpEyrcjc3Mh55hku79pF7uOPg6ritXQpwd262S+nqdgTAgkhRIUhZVoJ2AIDyZw1i5TvvqOwVSu0yckYnnwSwxNPoE1IcHQ8IYSo9KRMKxFLs2akrl9PxsyZ2Ly9cfvhB4K6d8fjiy9kK1UIIcqQlGllo9WSN3w4l3/+mfz770eTk4Pfc89hGDYMzaVLjk4nhBCVkpRpJWULCSF9yRLS33sPm68vbtu2Edy9O27ffOPoaEIIUelImVZmikJ+//5c3rYNU7duaDIyMIweje9zz6Hk5zs6nRBCVBpSplWArXp1jJ9+SsYbb6C6ueH5xRcE9u4t16UKIUQpkTKtKhSFvKFDSdmwAXN4OC6nThF43324r17t6GRCCFHhSZlWMZYmTUj97jvy/u//0OTn4z9pEr5TpshsNEIIcRukTKsg1dOTjHfeIf3tt+27fT/7jMCBA9FcvuzoaEIIUSFJmVZVikL+oEGkrl2LtXp19AcPEnTffbgcOeLoZEIIUeFImVZx5hYtSPnuOwratEGbnExg//64r1nj6FhCCFGhSJkKbEFBpK1eTe6QISgFBfiPG4fXO+/IqElCCFFMUqbCTq8nc/ZsMl95BVVR8Jk9235iktns6GRCCOH0pEzFVXKffJL0Dz6wn5j0+ecYnngCJSfH0bGEEMKpSZmKa5juu4/U1auxGgy4/fQTgQMGgIzrK4QQNyRlKq7LfOedpH7zDZawMFx+/RWX7t1lOjchhLgBKVNxQ9a6dUldvx5z06YoZ88S0K8f2thYR8cSQginoyuvFcXExLBs2TJsNhs9evSgX79+Vz2fmprKggULyM3NxWaz8eijjxIZGVle8cQN2AIDSf3yS6qNHIluzx4CBwwg7fPPsTRt6uhoQgjhNMply9Rms7F06VKmTp3K3Llz2bVrFxcvXrxqma+//poOHTowe/ZsJk6cyNKlS8sjmigG1dcXy6ZNmLp2RZuaSuDAgbgcPOjoWEII4TTKpUzPnj1LSEgI1apVQ6fT0bFjRw4cOHDVMoqikJeXB0BeXh7+/v7lEU0Ul6cnxuXLyb/vPjSZmQQMGSKFKoQQf1BUteyvzN+7dy8xMTGMGjUKgO3bt3PmzBlGjhxZtEx6ejqvvfYaubm5FBQU8NJLL1GvXr1r3mvr1q1s3boVgFmzZlF4iwO063Q6LBbLLb22LDl9LosF7YgRaFetQvX2xrJxI2r79o7P5WScNRc4bzbJVTK3k0uv15dyGlFux0xvZteuXXTr1o0+ffrw22+/MX/+fN566y00mqs3nqOiooiKiiq6n5qaekvrCwwMvOXXlqUKkWv2bPxMJjzWr0f7wAOkff455jvvdHwuJ+KsucB5s0mukrmdXDVq1CjlNKJcdvMaDAbS0tKK7qelpWEwGK5a5scff6RDhw4ANGzYELPZTHZ2dnnEEyWl05Hx7rvk9e2LJifHvss3OtrRqYQQwmHKpUzDw8NJSkri8uXLWCwWdu/eTevWra9aJjAwkOPHjwNw8eJFzGYzPj4+5RFP3Io/CjX/wQfRZGcT8OijuBw+7OhUQlRahYWQkKAlOtqF7793Y+9exdGRxF+Uy25erVbLiBEjmDlzJjabjbvvvpvQ0FBWrVpFeHg4rVu3ZujQoSxevJhNmzYBMGbMGBRF/rE4NZ2O9PnzQVVx37CBgMceI/Wrr7A0buzoZEJUKPn5CklJGpKStCQlaUlO1v7x/Z+PpaZqUNU/fyc+8oiVOXMcGFpcpdyOmUZGRl5z3eigQYOKvq9VqxYzZsworziitFwp1MJC3DdvJuDRR0ldswZr3bqOTiaE0zCbITFRy++/a7lwQUd8vJb4+CtftRiN2pu+h0ajUq2aleBgK8HBNiIinOaUF4ETnYAkKjAXF9IXLkQzdCiuu3YRMHgwqWvXYpOTHEQVYrPZCzM2VsfZs/ZbbKy9MBMTtVitN97TpterhIRYqV79ys32t/tWgoJs6P7yG9t+AlI5fDBRLFKmonS4uWH86CMCBg9Gf/gwAY88QtqaNdgCAhydTIhSZbHA+fM6TpzQcfasC7Gx2qKvJtP1T0NRFJUaNSzUrm3943b190FBNjQyuGuFJmUqSo3q5UXaJ5/YR0g6eRLDkCGkrV6NKieSiQrKaFQ4flxh715PTpxw4eRJHb/95oLJdP2tzOBgK+HhFurXtxAebr+FhVmoWdOKq2s5hxflSspUlCrV35+0zz8nsH9/9MeOYXjiCdI++wzc3BwdTYh/lJWlcOSICzExemJiXDhyRE9S0pVjmb5XLRsaaqFxYzMNG9qL80p5+viU+Rg4wklJmYpSZwsOJm3lSgL79cN17178J0wg/f33kf1YwllYLPDrry5ER7tw+LC9PGNjXa5Zzt3dRvPm0KBBPk2amGnSxEKjRmYpTXENKVNRJqyhofZdvgMG4L5xI9Zq1ch65RWQy52EA5hMEBOjZ98+++3gQT25uVf/cafXqzRtaiYiopCICDMREWbq1bMQHBxIamqmg5KLikLKVJQZS5MmGJcuJWDIELyWLsVaowa5f4zPLERZKiyE6Gg927e7smePnpgYPYWFV/8hFxZmoU2bQiIiCmnVykzjxmZkyFpxq6RMRZkq7NSJ9HfewTBmDL4zZmALCSH/b3PZCnG7VBViY7Vs3+7G9u2u7N599Zanoqg0bmymXbtC2rUroF27QqpVszkwsahspExFmTP17UtmUhK+M2bgN3Ei1sBACjt3dnQsUcGZTLBrlytbtrjx00+uJCRc/eusYUMzXbsW0KlTAW3bFuLnJ8c5RdmRMhXlIvfpp9EmJeH14YcYRo4kdc0aLE2bOjqWqGCMRoVt29zYssWNn392JS/vz61Pg8FK164FdO1aQJcuBdSoIVueovxImYryoShkTZuG9tIl+zi+Q4eSsnEjturVHZ1MOLmUFA2bNrmxYYM7+/frsdn+PPbZrFkhvXoV0LOniWbNzHLCuHAYKVNRfjQa0ufNQ3PpEq7792MYPpy0NWtQPT0dnUw4mfR0he+/d2f9end27fqzQHU6lc6dTfTqZaJXrwJq1rQ6OKkQdlKmony5uZG+dCmBffqgP34cv7FjSf/wQ9DefKBvUbmZTLB+vRtff+3B9u2umM32AnVxUene3cSDD+bTs6dJrvEUTknKVJQ7m8FA2ooVBPXti/uWLVhfe42sadMcHUs4gKrCsWMurFzpwfr1LmRkGADQalW6djXRt28+995rkpOHhNOTMhUOYa1fH+MHHxDwyCN4LVmCpW5d8oYOdXQsUU6MRg1ff+3OqlUenDz558hDLVoU8vDDefTpYyIwUE4gEhWHlKlwmMKOHcmYPRv/yZPx/e9/sdauTUG3bo6OJcrQkSMuLFvmyTffuFNQYN+NazBYGTAgn1GjXKleXeYUExWTlKlwqPxBg9CdP4/3/Pn4P/00qevXY2nUyNGxRCkqKICNG91ZtsyTw4ftQwwpiv046COP5BEVZUKvl/k5RcUmZSocLnvKFHRxcbhv2IBh6FBSN27EFhzs6FjiNqWlaVi2zJOPP/YgLc1+gpmfn41Bg/IYOjSXsDA5E1dUHlKmwvE0GtLnzkWbkIA+Oto+bdtXX6G6uzs6mbgFv/+uZfFiL1at8iia97NpUzNPPJFLv375uLvLyUSi8pEyFc7B3R3jsmUE9u6NPiYGvwkTSF+0SKZtq0COHXNh4UIvNm50K7ouNCrKxOjRObRrVygTBolKTX5TCadhCwzEuGIFNm9v3Ddtwnv2bEdHEsVw7JgLw4YZuPfeIL75xh2NBgYOzOPHHy+zYoWR9u2lSEXlJ1umwqlY7riD9MWLMTz+ON7z52OpV4/8hx92dCxxHceP63j7bW82b7bvjnd3t/HYY3n861851Kwpl7WIqkXKVDidgrvuInPGDPymTsVvyhSstWtT2L69o2OJP5w6peOtt7z59lt7ibq52Rg+PI/Ro3Pk2lBRZUmZCqeUN2wYuthYvJYuxTByJCkbN2KtW9fRsaq0pCQNc+Z4s3q1BzabgpubymOP5fLMMzkEB0uJiqpNylQ4raxp09DFxeG2bZt9lpkNG1D9/Bwdq8rJyVFYuNCLxYs9MZk06HQqQ4fmMm5cNiEhUqJCgJyAJJyZVkv6woWYGzdGd+4chqeeArPZ0amqDIsFPv7Yg06dgnnnHW9MJg3335/PTz9dZubMTClSIf5CylQ4NdXLC+OKFViDgnDdtQvfqVPto6OLMrVvn5577w3iP//xIzVVy513FrJuXQoffJBOvXoy2IIQfydlKpyetWZNjMuWobq54fn553guXuzoSJXWpUsaxo3zY8CAQE6edCE01MLixUbWr0+lTRvZKyDEjUiZigrB3KoV6e+8A4DPa6/h9v33Dk5UuZjNsGSJJ127BrNmjQeuriqTJmXz00+X6d3bJNeJimLTarVERETQrFkzBg4cSF5eXqm+//Llyxk7duw/LvPzzz+ze/fuovuLFi3i448/LtUcfydlKioMU+/eZD3/PIqq4jd2LC7Hjjk6UqUQE+PCffcF8corvuTkaIiKMvHjj5d57rlsZERHUVLu7u7ExMRw/Phx9Ho9ixYtKvcMfy/TUaNGMbSMp3gstzKNiYlhwoQJjBs3jnXr1l13md27dzNp0iQmT57MO39shQjxVznjxpE3cCCa/HwMw4dDQoKjI1VYeXkKr7ziQ58+9l26tWtbWL48jRUrjDIIvSgVXbp04ezZswC8/fbbNGvWjGbNmjFv3jwA4uLiaNSoEUOGDKFx48Y89NBDRVuyYWFhpP4xjdDBgwfpdp3pGTds2EC7du1o1aoVUVFRXLp0ibi4OBYtWsTcuXOJiIhgx44dTJ8+nTlz5gD2Lmrfvj0tWrSgf//+pKenA9CtWzeef/552rZtS8OGDdmxY0eJPmuxytRqtfLqq6+ycuXKEr35FTabjaVLlzJ16lTmzp3Lrl27uHjx4lXLJCUlsW7dOmbMmMHbb7/N8OHDb2ldopJTFDJmz6agfXu0ycnoBgxAyc11dKoKZ8cOPT16BLFkiRcAo0bl8OOPKfTsWeDgZKKysFgsfPfddzRv3pxDhw6xbNky9u3bx969e/nggw84fPgwAKdPn2bMmDGcPHkSHx8fFi5cWOx1dO7cmb1793L48GEGDx7M7NmzCQsLY9SoUUyaNImYmBi6dOly1WuGDh3Km2++ydGjR2nevDmvvPLKVZn379/PvHnzrnq8OIpVplqtloSEBDIyMkr05lecPXuWkJAQqlWrhk6no2PHjhw4cOCqZbZt28Y999yDl5f9f25fX99bWpeoAvR6jB98gCUsDE1MDH7jxoFVtqSKIzNT4emntQweHEh8vI7Gjc1s3JjKSy9lyWwuolTk5+cTERFB69atqV27NiNHjmTnzp30798fT09PvLy8GDBgQNGWX2hoKJ06dQLgscceY+fOncVe18WLF7nnnnto3rw5//vf//j111//cfnMzEwyMjK46667ABg2bBjbt28ven7AgAEA3HnnncTFxZXkYxd/N+9DDz3EgQMHOHHiBBaLpUQrMRqNBAQEFN0PCAjAaDRetUxiYiJJSUm89NJLvPjii8TExJRoHaJqUQ0G0lasQPXzw33zZnxef93RkZzezp16evQIZvlyLXq9yvPPZ/Hddym0bCln6YrSc+WYaUxMDPPnz0ev1//j8srfzm67cl+n02Gz2a9lNplM133tuHHjGDt2LMeOHWPx4sU3XK64XF1dAfsGZEl7rtgjIH344YcA12z6Kopyy7t//8pms5GUlMS0adMwGo1MmzaNOXPm4OnpedVyW7duZevWrQDMmjWLwMDAW1qfTqe75deWJclVAoGBqF99Bffei9eiRbi1bIltxAhHpwKc6+eVnw8vvaRl/nz7BN1t26p88IGZRo3cADfHhvsLZ/qZ/ZXkun1dunRh+PDhvPDCC6iqytq1a/nkk08AiI+PZ8+ePXTo0IHPP/+czp07A/ZjpocOHeK+++7j66+/vu77ZmZmUrNmTQBWrFhR9Li3tzdZWVnXLO/r64u/vz87duygS5cufPLJJ0VbqbfrtocTVItxAb3BYCAtLa3oflpaGgaD4ZplGjRogE6nIzg4mOrVq5OUlET9+vWvWi4qKoqoqKii+1cOUJdUYGDgLb+2LEmukgns0oW8WbPwe+45tOPGkeHvT+HfjpE4JJeT/LyOHXNh/Hg/fvtNi1Zrv9zllVfcyMhIxQniXcVZfmZ/Vxlz1ahRo5TT/LPIyEiGDx9O27ZtAXjyySdp1aoVcXFx3HHHHSxYsIARI0bQpEkTRo8eDcC0adMYOXIkL7300nVPPgKYPn06AwcOxN/fn+7du3P+/HkA+vTpw0MPPcT69euZP3/+Va9ZsWIFo0aNIi8vj3r16rFs2bJS+YyKWpw2BFJSUm74XFBQ0D++1mq1MmHCBF5++WUMBgP/+c9/GD9+PKGhoUXLxMTEsHPnTsaOHUtWVhbPP/88s2fPxtvb+x/fOzExsTjxr1EZ/wcpS86ey+e11/B6/31svr6kfvMNlr/9EeaoXI5is8GCBV7MmeONxaIQHm5m/vwMWrY0OzzbjUiukqlIZXojcXFx9O7dm+PHjzs6ym0r9pbplcJMTEwkOTmZyMjIYq9Eq9UyYsQIZs6cic1m4+677yY0NJRVq1YRHh5O69atadmyJUeOHGHSpEloNBoee+yxmxapEFdkTZ2K9vx53L//HsOwYaRu2IDtb3s/qoqUFA3jx/uxfbt9F+6IETlMnZotJxgJUYaKvWWak5PD3LlzOX78OIqi8O677zJ+/Hj69+/PoEGDyjrnDcmWafmoCLmUvDwCBgxAf+wYBe3akfbFF/DHCQWOzFWeduzQM26cPykpWgICrLzzTgZ333315S4V4b+lM6mMuZxly7QyKfbZvJ988gnHjx9Hp9OhqirBwcE0aNCAQ4cOlWU+IYpN9fDAuGwZ1pAQXPftw2/KlCozKL7FAv/7nzePPBJASoqWDh0K2LIl5ZoiFUKUjWKXaUxMDBEREfTs2bPosVq1anHp0qUyCSbErbBVr07aihXY3N3x+OorvN5919GRylxysoZBgwKYN89+WOTZZ7NYtSpNpkgTohwVu0wLCwuvuUwlOzsbnU7mFxfOxdKsGekLF6IqCj6zZ+P+1VeOjlRm9u+3T5W2d68rwcFWVq1KY/LkHLRaRycTomopdpnWqVOHQ4cOFY2z+PHHH3Po0CHCwsLKKpsQt6ygVy+ypk8HwO/ZZ3H95RfHBiplqgrLlnkwcKB9t27HjgX88EMKnToVOjqaEFVSsct08ODBWCwWzpw5A8CmTZtQFIWBAweWWTghbkfuk0+SM3o0isWC/5NP4nL0qKMjlYr8fJg40Y///tcPi0Xh6adz+OKLNAIDZbeuENczc+ZMmjZtSosWLYiIiGDfvn2lvo5i76Nt1KgRb7zxBlu2bCE1NZWgoCCioqKoU6dOqYcSorRkTZ2K5tIlPNaswfD446SuX4+1Au9NuXhRy5NP+nPsmB43Nxtvv51B3763N4SaEJXZnj172LhxI9HR0bi6upKamkphYenvwSnRAc/atWvz+OOPYzQaMRgMReMYCuG0NBoy3noLTWoqbtu3EzBkCKnr12OrIMOw/dW+fXpGjvQnPV1LnToWPvzQSJMmJRs/VIiqJikpicDAwKK+KqshGIu9mzcnJ4e3336boUOHMnHiRIYOHcrbb79NTk5OmQQTotTo9aR/8AGFzZqhi4vDMHRohZu2bfVqdwYNCiA9XUu3biY2bUqRIhUVj6KUze0f9OrViwsXLtCwYUPGjBnDL2V0/kSxy3TRokXX7Gfet2+fQ2ZRF6KkVC8vjJ98gqV2bfRHjuD/9NNgdv7ZUmw2eOMNbyZN8sdsVhg5MocVK4z4+1eN62eFuF1eXl4cOnSIJUuWEBQUxKBBg1i+fHmpr6fYu3mPHTtGcHAwzz33HDVr1uTixYvMmTOHY8eOlXooIcqCLTiYtE8/JbBfP9x++gm/554jY+5c0BT7b8pylZenMGGCH99+645WqzJjRibDhuU5OpYQt85Bg6hotVq6detGt27daN68OStWrGD48OGluo5i/xapVq0aTZs2pU6dOuh0OsLCwmjatCkhISGlGkiIsmQND8f4l0EdfKZNc8pRkpKTNfzf/wXw7bfu+PjY+PRToxSpELfg9OnTRVehgH0AorI4cfYft0xPnDhR9H3Xrl356quvCAsLo0aNGiQkJLBv3z65NEZUOObISNKXLsUwfDheH32E6u1N9pQpjo5V5NQpHY89FkBSkv1Eo48/NlK/vhwfFeJW5OTkMG7cODIyMtDpdNSvX58lS5aU+nr+sUz/PhE4cM3cb5988gkPPPBA6aYSoowV3HUX6QsX4v/003i/8w42Hx9yR41ydCz27tUzYoSBzEwNbdoU8NFH6RgMcv2oELfqzjvvZPfu3WW+nn8s04oyi7sQt8J0331kvPUW/hMn4jtjBqq3N3lDhjgsz6ZNbowb509BgcJ99+Uzf3467u4OiyOEKIF/LNMFCxaUVw4hHCJ/4ECU3Fz8XnwR3+efx+blhalv33LPsWyZBy+95IuqKgwblsuMGZkyvq4QFUiJR6nPzs6moODa+RGFqKjyhg9Hk5WFz5tv4j9+PEYPDwr+MjtSWVJVmDXLm/fes8/4MmVKFuPH59zs0jkhhJMpdpkePXqU999/H6PReNXjiqKwcuXKUg8mRHnKGTcOJTsb74ULMTz1FMalSyno3r1M12mxwL//7cfq1R5otSr/+18Ggwbll+k6hRBlo9iXxnzwwQfXFCmA6oSXFQhRYopC9tSp5IwYgVJYiOHJJ3H96acyW11BAYwa5c/q1R64u9tYtswoRSpEBVbsLdPs7GxatmzJ5MmTcXNzK8tMQjiGopD16qugqngtW4Zh5Ej7Furdd5fqavLyFEaO9Gfu7YZpAAAgAElEQVT7djd8fW18/HEarVs7/2hMQogbK/aW6b333ktqaipGo1G2RkXlpShkzZhB7rBhKAUFGEaOxPXnn0vt7TMzFR591MD27W4EBlr58stUKVIhypiXl1eZr6PYW6bt27dn8+bNTJo06arH5ZipqHQUhcyZM0FV8fz4YwwjRmD86CMKunW7rbdNS9Pw6KMGjh/XU726lZUrU6lf31o6mYUQDlXsLdN3332XvLxrhzOTrVRRKf1RqLmPP27fQh0xAtdt22757ZKSNAwYEMDx43rCwiysWydFKkRlUuwt09TUVMLDw3nsscfw8PAoy0xCOAeNhszXXwfA85NPMIwYQfr8+ZgefLBEbxMXp2Xw4AAuXNDRqJGZL75IIzhYRjUSVY/yStlc86VOc/xGXbHLNCoqirNnz9KwYUN0uhJfnipExaTRkPnGG6geHngtXoz/M8+QkZdH/uDBxXp5bKyWhx8OJDlZS6tWhXzySZpMnyZEJVSiKdguXLjAyJEjCQ4ORvOXaavefPPNMgknhFNQFLJeegmbtzc+c+bg/+yzaLKzyf3Xv/7xZWfP6hg4MIDLl7W0b1/AihVGvLykSEXV5QxbkGWl2GUaHx8PgMlkKvpeiCpDUciZNAnV2xvfadPwnT4dJSeHnIkTud5wRSdPwkMPBZCSoqVjR3uRenhU3l8kQlR1xS7T0aNHl2UOISqE3CefxOblhd+//43PnDlo0tPJmjaNvw6ke+qUjkcecSElRaFz5wKWLzfi7i5FKkRlVuwy7XablwUIUVnkDx6M6umJ//jxeC1dijYpifR33wV3d06e1PHwwwEYjQp33WVi6VKjzPwihIPl5OSU+TqKXaYLFy687uOKoshWq6hyTH36kBYQgGHkSNy//RZNSgo7nv+ch/9Vl/R0Lb162Xj/fSMyWJgQVUOxy/SXX3654XNSpqIqKuzYkdR16zA89hgnDpgZ9HAA6TYt3bub+PJLDeXwx7AQwkkUu0wfeuihou9tNhvx8fEcPHiQu0t53FIhKhLLHXew7c2tPDKsBuk2X3rrv+f9sTbc3KKkTIWoQopdpgMHDrzmscWLF193JpnriYmJYdmyZdhsNnr06EG/fv2uu9zevXt5++23eeONNwgPDy9uPCEc4vhxHYPG3kGmTUOfgB18lfYgLo9osL7/Ptxzj6PjCSHKSbGHE0xNTb3qFh8fT2JiIqdPn77pa202G0uXLmXq1KnMnTuXXbt2cfHixWuWy8/P57vvvqNBgwYl+xRCOMDJkzoGDw4gM1PDvffm8+6e2piHPoJSUIBuxAh8XnsNrDJkoBBVQbG3TJ955pnrPl67du2bvvbs2bOEhIRQrVo1ADp27MiBAweoVavWVcutWrWKvn378s033xQ3lhAOcfasvUjT07X06GHi/ffT0etdyHzjDcyNGuH78st4vf8+utOnSV+wANXHx9GRhRBlqNhbpn+n1+u54447bliyf2U0GgkICCi6HxAQcM3u4XPnzpGamkpkZOStRhKiXJw7p+XhhwNITdVy110mliwxotf/+XzesGFYvv0Wq78/bj/+SGDv3uiKsQdHCFE2tFotERERNGvWjD59+pCRkVHq6yj2lumqVatKfeVX2Gw2Pv74Y8aMGXPTZbdu3crWrVsBmDVrFoGBgbe0Tp1Od8uvLUuSq2TKO9f58/DIIy5cuqTQrZuNtWs1eHhcu35tjx5Ydu9GeeghXH79laDevbG+9x62IUPKLeuNyH/LkpFcFZ+7uzsxMTEADBs2jAULFvDiiy+W6jpuWqaDBg36x+eLM5+pwWAgLS2t6H5aWhoGg6Hovslk4sKFC7zyyisAZGRkMHv2bKZMmXLNSUhRUVFERUUV3U9NTb3ZR7iuwMDAW35tWZJcJVOeuRIStPzf/wVw8aJCmzYFLFliJC9P5TozE9pz+figrFuH7wsv4PH11+hGjCB361YyX30VR47kIP8tS6Yy5qpRo0Ypp6k4OnTowNGjR0v9fW95N+8VxZnPNDw8nKSkJC5fvozFYmH37t20bt266HkPDw+WLl3KggULWLBgAQ0aNLhukQrhKMnJGh5+2D6Nmn32FyOenjf/t696eJDxzjtk/O9/qK6ueH7+OUF9+6I9d64cUgvhXBSlbG7FZbVa2bZtGw+WcBrF4rjplunfZ4QxmUx8//337NmzB4CwsLCbrkSr1TJixAhmzpyJzWbj7rvvJjQ0lFWrVhEeHn5VsQrhbFJSNAwaFEBcnI7mzQv57LM0vL1LMNauopD36KMUtmiB4emn7bt9e/Uia/p08oYMKdlvAyFEieXn5xMREUFCQgKNGzemZ8+epb4ORS3OpiVQWFjI5s2b+eabb8jKyqJ27do89NBDtGvXrtRDlURiYuItva4y7ropS1U1V1qahoEDAzh92oXGjc2sXp2KwXDz/2VulEvJysJ36lQ81q4FwNSzJxlz5mArx2NfVfW/5a2qjLmq2m5eLy8vcnJyyMvL45577mHgwIGMHz++VNdx0y1Ts9nMli1b+Oabb8jIyKBWrVqMGDGCDh06lGoQIZxNerrC4MH2Im3Y0MzKlWnFKtJ/ovr4kPHee5h69sTvP//B7YcfCOrenYw5cyjo1auUkgshrsfDw4N3332Xfv36MWbMGHS6Yp+De1M3PWY6duxYPv74YzIzM+nQoQODBg3CxcWFgwcPFt2EqGyyshSGDAngxAkX6ta1sHJlGoGBtlJ7f1Pfvlz+4QcKOnVCm5ZGwBNP4DduHJq/nKgnhCh9rVq1okWLFnzxxRel+r43reUr1+OoqsqePXuKjpVeUZyzeYWoSHJy7EV65IieOnUsrF6dSrVqpVekV9hq1iRt5Uo8P/wQ7zffxGPNGlx/+omsV14hf8AAOZYqRCn5+xRsGzZsKPV13LRM5TomUZXk5SkMHWogOlpPzZoWVq9Oo0aN0i/SIhoNuU89halXL/yefx7XnTvxHz8e9zVryHzjDazFGGFMCOF4Ny3TBQsWlEcOIRwuPx+GDzewb58rISFWvvwyjVq1ymdsXWtYGGkrV+K+ejW+r76K288/43r33eSMHk3OM8+gygzjQji1277OVIjKoKAA/vUvA7t2uRIcbGX16lTq1CnnQeoVhfxBg7j888/k9e+PYjLhPXcuQV274rZhAxTvxHshhANImYoqz2yG0aP9+eknNwwGK6tWpREe7rjZXmxBQWS89x6pa9dibtoUXWIihlGjCBg4EN3x4w7LJYS4MSlTUaVZLDB2rD+bN7vj52dj5co0Gja0ODoWAIVt25Ly3XdkzJqF1d8f1z17CL7nHvyeeQZtXJyj4wkh/kLKVFRZNhtMnuzHxo3ueHvb+PzzNJo2dY4iLaLVkvf441zesYOcf/0LVa/HY906gu+6C98XX0Rz+bKjEwohkDIVVZTNBs8/78vXX3vg4WHjk0/SaNnS7OhYN6T6+5M1fTqXd+wgb+BAsFrxXL6c4E6d8H79dTROOEKPEFWJlKmoclQVXn7Zh88/98TNTWXFCiNt2jhvkf6VtVYtMubNI2XrVvJ79UKTl4f3ggUEt2uHz8svo7nF4TWFELdHylRUKaoKr73mw7JlXuj1Kh99ZKRjx0JHxyoxS6NGpC9bRsrGjfZSNZnwWrqUap064TtlihxTFaKcSZmKKmXOHG8WLfJCp1NZssTIXXcVODrSbTG3akX6smVc3rKF/D59wGzG87PPCO7cGf8RI9Dv3i2X1AhRDqRMRZXx7rtezJvnjVarsnBhOj17Vuwi/StL06akL1pEys8/k/fww+DigvvmzQQOHEhQz564r1oFJpOjYwpRaUmZiiph8WJP3nzTB0VReeedDB54oHIWi6V+fTLmzuXS/v1kPfss1qAgXE6exH/yZKq1aYPPK6/AyZOOjilEpSNlKiq95cs9ePVVXwDeeiuD/v3zHZyo7NmCgsiZPJlL+/aRPm8ehc2aoTUa8VqyBH1EBAH9+uG+ejVKfuX/WQhRHqRMRaX2xRcevPiiHwCvv57BoEFVrDxcXckfOJDU778nZeNGcocMQfXywvXAAfwnTaJaq1b4Pvcc+p07weq4UZ+EqOikTEWl9dVX7vz73/Yt0mnTMhk2LM/BiRxIUTC3akXm7NmYf/+djDlzKGzVCk12Np5ffEHgoEH23cDTp+MSEyMnLQlRQlKmolL66it3Jk70Q1UVXnghi6eeynV0JOfh5UXeI4+QunEjl3/6iewJE7DUqYP20iW8PviAoAceILhzZ7xnzsTl4EH7CBdCiH8kZSoqna+//rNIp0zJYty4nJu/qIqyNGxI9pQpXN61i5QNG8gZORJrUBC6uDi8Fy4kqG9fqkVG4jtlCq7btskZwULcwE3nMxWiIlmz5s8i/fe/s5gwQYq0WBQFc2Qk5shIsl5+Gf3+/bht3ozb5s3oLlzA87PP8PzsM2yenhR06kTBXXdR0LUr1rp1QVEcnV4Ih5MyFZXGmjXuTJjgh81mL9KJE6VIb4lOR2HHjhR27EjW9OnoTpwoKlb98eO4b9mC+5YtAFhCQyno2pWCbt0o6NAB1d/fweGFcAwpU1EprF37Z5E+95wUaalRFCxNm5LTtCk5kyejSUzEdft23H75Bdft29FduIDuj61WAHOjRhS2a0dBu3YUtm2LrXp1B38AIcqHlKmo8NaudWf8+D+LdNIkKdKyYqtRg/zBg8kfPBisVlyOH8f1j2LVR0fjcuoULqdO4bliBQCW2rUpbNuWwshIzBERmBs1AldXB38KIUqflKmo0Fat0hQV6bPPSpGWK60Wc8uWmFu2JGf8eDCZ0B89in7fPvvt4EF08fHo4uPx+OorAFS9HnPjxphbtqQwIgJzixZYGjQAnfwqEhWb/AsWFdaqVe48+6wWVVWYPDmbyZOlSB3Kzc2+Fdq2LYwbB1YrupMncd2/H5eYGFyOHEEXG4v+yBH0R47g+fHHAKiurpgbNEAbEYFn3bpYGjfG3LgxtqAgOblJVBhSpqJCWr78z5GN5GQjJ6XVYmnWDEuzZkUPKdnZuBw9isvRo+ivFOyFC+iPH4fjx/H9y8utBgOWRo2w1K+PJTwcS716WMLDsdaqBVpt+X8eIf6BlKmocBYt8mTGDPuv3dmzLQwZIkVaUaje3hR26kRhp05cGUZDycrC5fRp/C5coODgQXQnT+Jy6hRaoxHt7t247t599Xvo9VjCwv4s1zp1sIaGYqlVC2vNmnJMVjiElKmoMFQV5s3zYs4cHwDeeCODCRM8SE11cDBxW1QfHwrbtMF2331kDhjwx4Mq2sREdKdOoTt3Dl1sbNFXbXIyLr/9hstvv133/awhIVhr1cISGoq1Vi2soaFYa9SwP16tmv3yHdl9LEqZlKmoEFQVXn/dm4ULvdFoVN56K4OHH84HPBwdTZQFRcFasybWmjUp6NHj6qdyc9GeP/9nwcbHo71wAe3Fi2gTE9EmJ6NNTkZ/8OB131p1dcVarRrWatWwXfkaEmIv26AgbAEB9pvBAC4u5fFpRSUgZSqcns0GL73ky/Llnuh0KvPnp/PggzKsXVWlenpecyy2iMViL9MLF4oKVhcfjzYpCc2lS2gvXUKTlVV0lvHN2Hx9UYKCCPT1xXqlZP8oWltAADZfX1Q/P2w+Pth8fbH5+oKbm2z5VkHlVqYxMTEsW7YMm81Gjx496Nev31XPb9y4kW3btqHVavHx8WH06NEEBQWVVzzhpAoKYOJEf775xh29XmXxYiO9ehU4OpZwVjqdfddurVrQocN1F1Hy8tAkJ6O9Uq5/bMlqk5PRpKaiMRrtX9PT0WRmQmYm+hJEUPV6e7H6+KD+UbA2X19UHx/7Y15e2Ly8UD08UL287PevfO/pWXQfV1cp5QqkXMrUZrOxdOlS/vvf/xIQEMB//vMfWrduTa1atYqWCQsLY9asWbi6urJlyxY+/fRTJk2aVB7xhJPKyVF48kkDO3a44uVlY+lSI507Fzo6lqjgVA8PrPXqYa1X758XtNlQMjIIsNnIio1Fk5b2581otH/NykKTkYGSlYUmMxNNZiZKYSHalBS0KSm3l1OnQ/X0xObpaS9ZT09Ud/eim+buu+GRR25rHaL0lEuZnj17lpCQEKpVqwZAx44dOXDgwFVl2uwvu2waNGjAjh07yiOacFKpqRqGDjVw5IieoCArn36aRrNmFkfHElWJRoNqMEBgIIWBgcV/XX6+vWQzM1H+KFhNVlbR90peHpqcHJTc3KKbJicHJS8P5Y/HNbm5KIWFRa+5HquXl5SpEymXMjUajQQEBBTdDwgI4MyZMzdc/scffyQiIqI8ogknFB+v5dFHAzh/XkedOhY+/zyNsDCro2MJUTzu7tjc3bH9sfFwywoL7cV6pWTz8uy3/HyUvDy8GzUqnbyiVDjdCUjbt2/n3LlzTJ8+/brPb926la1btwIwa9YsAkvyF+Nf6HS6W35tWarquWJiFPr315GcrBARYeObb2xUq3bjmUiq+s/rVjhrNslVMlqdjkCL7K1xFuVSpgaDgbS0tKL7aWlpGAyGa5Y7evQoa9euZfr06bjc4JT0qKgooqKiiu6n3uJFhoGBgbf82rJUlXNt2eLKM8/4k5en0LFjAR99ZESrVf/xOtKq/PO6Vc6aTXKVzO3kqlGjRimnEZryWEl4eDhJSUlcvnwZi8XC7t27ad269VXLnD9/ng8++IApU6bg6+t7g3cSldXSpZ6MHGkgL0/DgAF5fPppGt7eqqNjCSFEsZTLlqlWq2XEiBHMnDkTm83G3XffTWhoKKtWrSI8PJzWrVvz6aefYjKZePvttwH7X13PP/98ecQTDmSxwPTpPixb5gVQNBepXBEghKhIyu2YaWRkJJGRkVc9NmjQoKLvX3rppfKKIpxETo7C6NH+/PijG3q9yty5GfTrl+/oWEIIUWJOdwKSqBpiY7WMHGngzBkXDAYrH32UTps2cg2pEKJiKpdjpkL81datrvTuHcSZMy40bGhmw4ZUKVIhRIUmW6ai3Nhs8M47Xrz1ljeqqnD//fnMnZuBl1fZnWiUb8knOTeZ5LxkknOTSTelk1WYRbY5m+xC+81is2BVrVhVKzbVhqqquGpdcdO54aZ1w03nhqeLJwY3A/5u/gS4BWBwMxDkHkSIZwguGhkMXYiqTspUlIvsbIWJE/34/nt3FEXl+eezGDeudE40sqk2zmWe47f034jNjOVsxlliM2OJy4wjvSD99lfwDzSKhuqe1Qn1CqWWdy1CvUNp4NeAdrZ2+Kv+uGplbk0hqgIpU1Hmjh/X8fTTBuLidPj42HjvvXR69Lj1weoTchI4kHyAI6lHOJlxksPJh8kxX3+CcJ2io5pnNUI8QgjxDCHQPRBvvTc+Lj72r3ofXLQuaBUtGkWDVtECUGAtwGQ1YbLYbznmHIwmI2mmNNJN6aSZ0ricd5lLeZdIyEkgIScBkq9et1bREuYTxh3+d9AkoAktg1oSERSBwe3aa6yFEBWblKkoM6oKK1Z48MorvhQWKjRubGbJEiP16pVsaMCEnAR2J+5mT9Ie9iTtIT772qmzqntWp7GhMfX96hPuG064Xzj1fOsR5B6ERim7UwMKrYUk5iZyIfsCF7MvEpcdx5n0M8RmxXIu4xyxmbHEZsbybdy3Ra8J9QqlRVALIoIiaFOtDS2DWqLXlmReEiGEs5EyFWUiM1Phuef8+PZbdwAefzyXadMycXe/+WutNiuHLh/ih99/YGv8Vn7L+O2q5330PrQNaUtEUARdwrsQpg8j0N0xw73ptXrCfMII8wm76vHAwEAuJF/gXOY5Tqef5ljqMY6kHOFo6lEu5FzgQs4FNp3fBICb1o3I4Eg6VO9Au+rtiAyOxF1XjB+UEMJpSJmKUrd/v57x4/24cEGHt7eN2bMzbjqZd74ln23x29j8+2Z+vPAjGQUZRc95uXjRvnp7OlbvSMcaHWliaIJWY98d66xDvQG469xpGtCUpgFNGVB/AGD/Q+FMxhmOpBwh+nI0+5P381vGb+xO2s3upN0AuGhcaBXUii41u9ClZhcigiPkJCchnJyUqSg1JhPMmePDokWeqKpCixaFvP9++g1nfLHarOxO2s2as2v49vy3Vx33DPMJI6p2FD1r96RtSNtKsxtUq9HSyNCIRoZGDLrDPmhJWn4a+5L3sTd5L/uS9vFr2q/sv7Sf/Zf281b0W3i5eNGhege61uxKl5pdqO9XH0WGiBLCqUiZilJx/LiOCRP8OXXKBY1GZfz4bCZOzEb/tw5UVZVfjb+y5swa1seuJznvz7N2Wga2pE+9PvSs05Nw3/AqUxgB7gHcX/d+7q97PwCZBZnsTdrLjoQdbE/YTmxmLD/E/8AP8T8AEOIZQpcaXehaqyuda3Qm2CPYkfGFEEiZittUWAgLF3oxb543ZrNC3boW3nknnTvvNF+1XEJOAmvPrmXN2TWcTj9d9Hht79oMqD+A/vX7U9+vfnnHd0q+rr7cE3YP94TdA9h/djsTdrI9YTs7EnaQnJvMl2e+5MszXwLQ2NCYrjW70rVmV9pVbyfHW4VwAClTccsOHHDhhRf8OHXKfjzviSdyePHFbNzd7YMwZBRksOn8JtacWcPe5L1Fr/N39efB8AcZUH8AdwbfWWW2QG9VTa+aDLpjEIPuGIRNtXHSeJIdCTvYkbCDvUl7OWk8yUnjSRYfW4yr1pU21dpwV6276FqzK00CmpTp2cxCCDspU1FimZkKb7zhw6efeqCqCmFhFt58M4POnQspsBbw7fltrD27lq3xWym02YcJdNO60atOL/rX70+3Wt0qzTHQ8qZRNEUnNY1qMQqTxcTBSwfZnrCd7QnbOZZ6jJ2JO9mZuJOZzCTALYAuNbsUHW+t4SXzWApRFqRMRbGpKqxb586MGT5cuqRFp1MZMyabseMyOZq5jyk71rLx3EYyCzMBUFDoXKMzAxoM4P6w+/HWezv4E1Q+bjo3OtfsTOeanZnKVNLy09iZuJPtF7fzS8IvJOUmsS52Heti1wHQwK8B99S/hzYBbehQvQOeLp4O/gRCVA5SpqJYDhxQGD8+kOho+xZl69aFjH7xMIf5lG7r1pCYm1i07JVLQfqG96W6Z3VHRa6SAtwD6Bvel77hfVFVldjMWH65+AvbE7azJ2kPZzLOcObgGcB+CU7raq3tx1trdaV5QPOiS46EECWjqKpadqOMl4PExMSbL3Qdznp9orPlSkrS8MYbPnz9tQcAAYFmOj62kdi6L3Mi/XjRcjW9atK/fn8GhA/gDsMd5ZbP2X5eVzhjrkJrIdGXozloPMj3Z77nSOoRbKqt6Hk/Vz861+hcdLy1lnetcs3njD8zqJy5atSQ3f2lTbZMxXUZjQrvv+/FRx95YjJp0LlYCenxBRcjnmGDWxakg6/el971ejOg/gDahrSVE12cnF6rp3319vRu3puxTceSbkpnV+Iu+/HWi9u5kHOBjec3svH8RgDq+tQtKtaONTrKbnoh/oGUqbhKdrbChx96smixJznZ9l1+SpO1WKKe46LhHHqNnqja9zOg/gC61+4us6JUYP5u/vSu15ve9XqjqipxWXFFxborcRfns85z/sR5lp9YjkbR0DygOe2qt6N9SHvahLSRAfuF+AspUwHYz9BdvsKdhYvcycn8oyDrbYHu/4VaB+lW5y4eqD2a++vej5+rn2PDilKnKAp1fetS17cuw5oMw2KzcDjlMDsu7uCXhF+IuRzDkdQjHEk9wpJjSwBo5N+IdtXb0S6kHe2rt6eaRzUHfwohHEfKtIpLSFSZ8U42339dB3P+Hxf7h+6EHi/SonUG/cL78WD4+zSv09wpjxuJsqHT6GhTrQ1tqrVh8p2TyTXncujyIfYl7WNf8j6iL0dzKv0Up9JPseLECsA+BGS7kHbcWe1OIoMjaejXUE5oElWGlGkVVGAtYOX243y41JNz2zuBtab9ibrbCOm1gkfuD6Jf/RkyIpEo4uniWTTKEtj/DR1JOcLepL3sS97HgUsHiMuKIy4rjlW/rSp6TURQBJHBkUU3R83uI0RZkzKtIrIKs/jx/E6Wf51J9Kb2WOP6/PGMDY8W3/HAY8d5olczWgS+ISMSiZty1brSNqQtbUPaAmCxWewD9Cfv53DKYaIvRXMh5wK7EnexK3FX0etqe9cmMjiSFoEtaB7YnKYBTfF19XXUxxCi1EiZVlKqqnIm4ww/XviRTYdOcviHlqiHRkKO/bpPjWsuzbsfYdJoDVGRLVCUlg5OLCoynUZHy6CWtAz689/R5bzLHL58mOjL0Ry6fIiYlBjis+OJz44vGkQCoI53HZoFNqN97fbUc69Hs8BmsgUrKhwp00rkymwjvyT8wtbTB0nY3wGODIX4V4uWCax9iaHDcnhqiAfe3mGOCysqvWCP4KsG7LfYLJxOP0305WiOpR7j17RfOWk8ye/Zv/N79u9Fk6WDfWacpoamNDI04g7/O7jDcAf1fevjpnNz1McR4h9JmVZgOYU57L+0n92Ju9mVuIujCefg7L3w60A4vQis9l88ejcz996Xy+OPWujQwYqiuAMVeqwOUQHpNLqicYWvMNvMnEk/w7G0Y8TmxnLgwgF+Nf5Kcm4yybnJbLuwrWhZjaIhzCeMRv6NuMNwBw39GtLI0IgwnzAZ61k4nJRpBZKcm1y0y+xA8gGOpBzBkucFv/WGE69A7D1gsZ+RqygqnTqbeOihfO6/34Snp5SncD4uGheaBDShSUCTohF9bKqNc5nnOGk8yW/pv3Eq/RSnjac5n3Wec5nnOJd5jm/jvi16D42iobZ3ber51qOuT13q+dYrutXwqiGDiYhyIWXqpPIt+Zw0nrSX56VDHLp8iIScBLApkBxhL87YtyC+I9hcil53552F3H9/Pn36mKhZ0+rATyDErdEoGur71b/mbHKTxURsZiyn009z2nja/jX9NBeyLxSdSfx3rlpXwnzCqOdbj8CR/5QAAA+CSURBVNretQn1DqWWVy1CvUMJ9Q6Vgf5FqZEydQJGk5Hjacf5NfVXzuacJToxmrOZZ+3jpqpAel2I74Eu7l6U2HswZ/858oxGo9K+YwEPPJDPvfeaCAmx3XhFQlRgbjq3a3YTg71k47PjOZd5jvOZf269ns86z6W8S0Wlez3+rv7U9q5NLe8/CtYrlBpeNQjxCCHEM0ROhBLFJmVaTlT1/9u795im7v+P488W5A4tbaGAClEu3+/mz6+X1UtQ51RmlmV/GKNuLt8YtywmP3ROjdkkmWY3MxNHWOIl++O3Eef2h26R+M9+2eKcmHmJiOIWppvg3OoXaGnLpRQovZzfHx39ioJfsXBaf7wfSdNz2tPy4nNOz7vnc07PUWjrbaOps+m+W1tv278nDGqh/Uk0f/436S3P4b9dSl9HqHj6/5okNzfAM8/088wzXhYv9qLXSxeumLiS4pMoySyhJLPkvud6BnrC3cNWtxWr28qdnjvh+w5vBx3eDq45rg373vGaeHLTc8lKyiInNYfclFxyUnMwp5jD98YkI/pEvfykbIKTYjqG+v39/KvnX6EPbY+VO+47/On+kz+6/6CpqwmPzzP0BUEtOEtIsD2P3rUCWp6i689peHsTUAD3X5Pp9UHmz/eycOEAy5Z5KS72I59bIf6ztIQ0ZppmMtM0877ngkqQ9r72UGF138HaEyq2LT0ttPWGDoDq8HZg7bZi7bY+8O/Ea+IxJhsxJhkxJZswJZuGHTYkGdAl6shIyJB9uf/PSDF9CIqi4Pa5sXls2Hrvv7V6WrG6rdh6bSO8AdA9hbSu5zB0L2KSczbetkLab+fi7ZvEAGC/a/LJk/3Mnz/A/PkDLFgwQHGxH6187oQYU1qNFnOKGXOKGYvZMuw0/f5+BhIHuH7nOm29bbR6WmnztGHrtdHmaaO9rx1nv5Puge7w+uBhaNCgS9ShT9SjT9SjS9ChT7pr+K/H9Yl6MhIzSE9IJ2NS6D49IZ14ray6Y41qc6ShoYHq6mqCwSArVqxg1apVQ573+XwcPHiQW7dukZ6ezrZt28jOzh6XLN6AlxZ3C03OplA3T3+oq8fV7woP3/1Ym6eN/kD/g980qEHbk4+pfz4Zntkkuv+G4pqGxzYZ2x9G+nsn0QP03POyvDw///iHj5kzfcya5WPp0nS0WjkHrhCxICk+iSmZU8gIZDxwOm/Ai7PPibPfiaPPgaPPMexwR38Hnd5O3D43nd5OOr2dj5QrOT6ZtU+s5cOFHz7S68XYU6WYBoNBPv30U95++22MRiMVFRVYLBamTPn3xYdPnz5NamoqBw4c4Ny5c3z55Zds3759XPL883//yfnW8w//gqCWZO80Mgf+i3Tv30jqm4bWnY/izsXbYaKzzYCjNRWfT4udoVuZg4zGACUlfoqL/ZSU+Cgu9vP3v/sxmYYeMGQypSPnkxfi8ZIYl0heWh55aQ930W1/0E/3QHe4oHZ6O+nydg0ZH7x1e7tx+9y4B/66+dz0+fsIKHK0fixRpZg2NTWRk5OD2Ry6RFNpaSl1dXVDiunly5dZu3YtAAsXLuSzzz5DUZRx2alvjJ+MyT+LVG8+yb6pJA1MIb4vB21/ForHSMCjZ6Ang77uVLo7knG1J9AX0ND3H97XZAqQnx8gP99Pfn6AggI/BQWhImo0ylG2QoiQeG08hiTDI10TVlEUPD4PmYZMfD2+cUgnHoUqxdTlcmE0GsPjRqORmzdvjjhNXFwcKSkpuN1uMjIe3L3yKBz/cxTHhURGswFoMAQwm4Pk5gYwm0PDOTmh4alTQ0VUTowghBhvGo2GtIQ0dEk6HD3SjRUrHru92KdOneLUqVMA7Nu3D5Np9L8DKyiI4/ffFYxGBZMJDAbuGh56n5WlkJsLiYmDr9b+dRsf8fHxj/Q/jTfJNTqxmgtiN5vkGp1YzTVRqVJMDQYDTqczPO50OjEYDMNOYzQaCQQC9Pb2kp6eft97lZWVUVZWFh5/lAtWV1YSPnXZw3C7Qzc1jCaXmiTX6MRqLojdbJJrdCLJlZf3cPt2xcNT5QcXhYWFtLa2Yrfb8fv9nD9/Hotl6KHoTz31FGfOnAHg4sWLzJgxQ34ELYQQ4rGgypZpXFwcr776Knv37iUYDLJs2TKmTp3KsWPHKCwsxGKxsHz5cg4ePMjrr79OWloa27ZtUyOaEEIIETHV9pnOnTuXuXPnDnnsxRdfDA8nJCSwY8cOteIIIYQQY0bOqyOEEEJESIqpEEIIESEppkIIIUSEpJgKIYQQEZJiKoQQQkRIoyiKnANPCCGEiMCE3TLdtWtXtCMMS3KNjuQavVjNJrlGJ1ZzTVQTtpgKIYQQY0WKqRBCCBGhuHfeeeedaIeIlunTp0c7wrAk1+hIrtGL1WySa3RiNddEJAcgCSGEEBGSbl4hhBAiQo/dxcHHQkNDA9XV1QSDQVasWMGqVauiksPhcHDo0CE6OzvRaDSUlZXx/PPPc/z4cb7//nsyMjIAWL9+/X0XCRhvmzdvJikpCa1WS1xcHPv27aOnp4eqqira29vJyspi+/btpKWlqZappaWFqqqq8LjdbmfdunV4PB7V2+vw4cNcuXIFnU5HZWUlwIjtoygK1dXVXL16lcTERMrLy8ete264XEePHqW+vp74+HjMZjPl5eWkpqZit9vZvn17+NqWxcXFbNq0SbVcD1rOa2pqOH36NFqtlldeeYXZs2erlquqqoqWlhYAent7SUlJYf/+/aq210jrhlhYxsQIlAkmEAgoW7ZsUdra2hSfz6fs3LlTsVqtUcnicrmU5uZmRVEUpbe3V9m6datitVqVY8eOKSdPnoxKpkHl5eVKV1fXkMeOHj2q1NTUKIqiKDU1NcrRo0ejEU1RlNB8fO211xS73R6V9mpsbFSam5uVHTt2hB8bqX3q6+uVvXv3KsFgUPn111+ViooKVXM1NDQofr8/nHEwl81mGzLdeBou10jzzWq1Kjt37lQGBgYUm82mbNmyRQkEAqrlutuRI0eUr776SlEUddtrpHVDLCxjYngTrpu3qamJnJwczGYz8fHxlJaWUldXF5UsmZmZ4W+PycnJTJ48GZfLFZUsD6Ouro6lS5cCsHTp0qi1G8DPP/9MTk4OWVlZUfn7Tz755H1b5SO1z+XLl3n66afRaDSUlJTg8Xjo6OhQLdesWbOIi4sDoKSkJCrL2HC5RlJXV0dpaSmTJk0iOzubnJwcmpqaVM+lKAoXLlxg0aJF4/K3H2SkdUMsLGNieBOum9flcmE0GsPjRqORmzdvRjFRiN1u5/fff6eoqIgbN27w7bffcvbsWaZPn86GDRtU7U4dtHfvXgCeffZZysrK6OrqIjMzEwC9Xk9XV5fqmQadO3duyEouFtprpPZxuVyYTKbwdEajEZfLFZ5WTadPn6a0tDQ8brfbefPNN0lOTuall17iiSeeUDXPcPPN5XJRXFwcnsZgMETlC8D169fR6XTk5uaGH4tGe929bngclrGJasIV01jU399PZWUlGzduJCUlhZUrV7JmzRoAjh07xueff055ebmqmd5//30MBgNdXV188MEH4f1EgzQaDRqNRtVMg/x+P/X19bz88ssAMdFe94pm+4zkxIkTxMXFsWTJEiC09XP48GHS09O5desW+/fvp7KykpSUFFXyxOJ8u9u9X9ii0V73rhvuFovL2EQ24bp5DQYDTqczPO50OjEYDFHL4/f7qaysZMmSJSxYsAAIfePUarVotVpWrFhBc3Oz6rkG20Sn0zFv3jyamprQ6XThrqOOjo7wgSNqu3r1KtOmTUOv1wOx0V7AiO1jMBhwOBzh6aKxzJ05c4b6+nq2bt0aXgFPmjSJ9PR0IPR7RbPZTGtrq2qZRppv935GXS6X6u0VCAS4dOnSkK14tdtruHVDLC9jE92EK6aFhYW0trZit9vx+/2cP38ei8USlSyKovDJJ58wefJkXnjhhfDjd+/ruHTpElOnTlU1V39/P319feHhn376ifz8fCwWC7W1tQDU1tYyb948VXMNuneLIdrtNWik9rFYLJw9exZFUfjtt99ISUlRtfutoaGBkydP8tZbb5GYmBh+vLu7m2AwCIDNZqO1tRWz2axarpHmm8Vi4fz58/h8Pux2O62trRQVFamWC0L75PPy8obsElKzvUZaN8TqMiYm6Ekbrly5wpEjRwgGgyxbtozVq1dHJceNGzfYs2cP+fn54a2F9evXc+7cOW7fvo1GoyErK4tNmzap+sGw2Wx89NFHQOgb+uLFi1m9ejVut5uqqiocDkdUfhoDoeJeXl7OwYMHw91eBw4cUL29Pv74Y3755Rfcbjc6nY5169Yxb968YdtHURQ+/fRTrl27RkJCAuXl5RQWFqqWq6amBr/fH55Xgz/puHjxIsePHycuLg6tVsvatWvH7YvlcLkaGxtHnG8nTpzghx9+QKvVsnHjRubMmaNaruXLl3Po0CGKi4tZuXJleFo122ukdUNxcXHUlzExvAlZTIUQQoixNOG6eYUQQoixJsVUCCGEiJAUUyGEECJCUkyFEEKICEkxFUIIISIkxVSIR3Do0CHWrVvH8ePHox1FCBED5HSCQgxj8+bNtLe3j/j81q1bSU1NpaSkRMVUQohYJcVUiGEsW7aMnp4eAL777jv8fj8LFiwInxGnqKiIxYsXRzOiECKGSDEVYhiDJ2CH0Hlt/X4/zz33HDNmzABC3by1tbWsWbMm3N379ddfM3v2bNLS0sKnx3vjjTc4efIkP/74I9nZ2WzevDl8aS2Hw8EXX3zB9evX8Xq9FBUVsWHDBvLz86PyPwshHp3sMxViDF27do2enh5MJhPNzc3s2rWL27dvU1BQgNVqpbq6GgCv18u7777LhQsXKCgowGKx0NjYyHvvvUd3d3eU/wshxGhJMRViDJnNZioqKli1ahUQKpq7d+9m06ZNANy+fRsInR/aZrORmZlJXl4eaWlpmEwmuru7uXjxYrTiCyEekXTzCjGG8vLy0Gg04ZPw6/V6UlJSSEpKAkLFFQgf3ORyufjmm2+GvEdbW5uKiYUQY0GKqRBjSKvVPnB8UFZWFhC6LuaHH34YvjKIx+NBrj0hxONHiqkQUTBnzhyys7O5desWu3fvJj8/H6fTSWNjIxUVFeEDnYQQjwfZZypEFCQlJbFnzx4WLVqEw+GgtraWlpYWlixZQl5eXrTjCSFGSa5nKoQQQkRItkyFEEKICEkxFUIIISIkxVQIIYSIkBRTIYQQIkJSTIUQQogISTEVQgghIiTFVAghhIiQFFMhhBAiQlJMhRBCiAj9H87PxnyPzuepAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sline = plt.plot(\"t\",\"S\",\"\",data=sir_out,color=\"red\",linewidth=2)\n", "iline = plt.plot(\"t\",\"I\",\"\",data=sir_out,color=\"green\",linewidth=2)\n", "rline = plt.plot(\"t\",\"R\",\"\",data=sir_out,color=\"blue\",linewidth=2)\n", "plt.xlabel(\"Time\",fontweight=\"bold\")\n", "plt.ylabel(\"Number\",fontweight=\"bold\")\n", "legend = plt.legend(title=\"Population\",loc=5,bbox_to_anchor=(1.25,0.5))\n", "frame = legend.get_frame()\n", "frame.set_facecolor(\"white\")\n", "frame.set_linewidth(0)" ] } ], "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.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }