SIR-case-study.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba5b1899",
   "metadata": {},
   "source": [
    "# S.I.R. Modell\n",
    "\n",
    "Diese Case-Study ist vom Blog-post von Bhaskar Krishnamachari [1] sowie der Fallstudie aus dem Buch \"Learning Scientific Python\" von Christian Hill [2] inspiriert. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "109bb164",
   "metadata": {},
   "source": [
    "Das S.I.R Model ist ein klassischer Ansatz in der Epidemiologie zur Modellierung der Ausbreitung von Krankheiten mit Immunität in einer Population. Jede Person wird dabei in eine von drei Gruppen eingeteilt:\n",
    "* Gruppe S (Susceptible): Nicht immune, gesunde Personen\n",
    "* Gruppe I (Infected): Infizierte Personen\n",
    "* Gruppe R (Recovered): Genesene (immunisierte) Personen\n",
    "\n",
    "Es werden folgende vereinfachenden Annahmen gemacht:\n",
    "* Jede Person wird nur einmal infiziert und wird danach gesund oder stirbt.\n",
    "* Infizierte sind sofort ansteckend.\n",
    "* Die Infektionsrate ist konstant.\n",
    "* Die verstorbenen Personen werden wie die Immunisierten zu R gezählt. \n",
    "\n",
    "Mathematisch entspricht das Modell einer gekoppelten Differentialgleichung. Wir werden aber gleich sehen, dass dies in unserer Simulation etwas ganz einfaches ist. Das Modell definiert dabei wie sich S, I und R über einen Zeitschritt ändern und ist wie folgt definiert:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{dS}{dt} &= -\\beta \\frac{SI}{N}  \\\\\n",
    "\\frac{dI}{dt} &= \\beta \\frac{SI}{N} - \\gamma I  \\\\\n",
    "\\frac{dR}{dt}  &= R_t + \\gamma I_t\n",
    "\\end{align}\n",
    "$$\n",
    "Dabei bezeichnet die berechnete Zahl $\\frac{dX}{dt}$ jeweils die Änderung der Grösse $X$ in einem Zeitschritt.  Die Zahl $\\beta$ bezeichnet die Ansteckungsrate (wieviele Ansteckungen ein Fall pro Zeiteinheit verursacht) und die Zahl $\\gamma$ die Rate mit der sich infizierte erholen. Die Zahl $N = S + I + R$ bezeichnet die totale Anzahl Personen in der Simulation.\n",
    "\n",
    "Wie erwähnt ist in unserer Computersimulation das Berechnen der Differentialgleichung sehr einfach. In unserer Simulation entspricht ein Zeitschritt einem Tag. Wir berechnen die zu erwartende Änderung pro Gruppe (S, I, R) während des Tages und addieren diese dann zur Situation am aktuellen Tag. So erhalten wir die Prognose für den Zustand am nächsten Tag. Wir berechnen also Folgendes:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "S_{t+1}  &= S_{t} + (-\\beta \\frac{S_{t}I_{t}}{N})  \\\\\n",
    "I_{t+1} &= I_{t} + \\beta \\frac{S_{t}I_{t}}{N} - \\gamma I_{t}  \\\\\n",
    "R_{t+1}  &= R_{t} + \\gamma I_{t}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Dabei bezeichen $S_{t}$, $I_{t}$ und $R_{t}$ jeweils die Situation am aktuellen Tag und $S_{t+1}, I_{t+1}, R_{t+1}$ die Situation am nächsten Tag. Sie sehen, alle Gleichungen bestehen aus einfachen Produkten, Brüchen, Additionen, Subtraktionen bzw. dem Umkehren des Vorzeichens.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5f7cad4",
   "metadata": {},
   "source": [
    "Wir werden dieses Modell nun in Python simulieren indem wir die Berechnungen implementieren und dann wiederholt in einer Schleife anwenden, um die Werte für eine bestimmte Anzahl Tage zu erhalten."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ecbe1dd",
   "metadata": {},
   "source": [
    "#### Bibliotheken\n",
    "\n",
    "Um die Verläufe zu plotten nutzen wir wieder pyplot, welches wir mit der nächsten Zeile importieren. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db75f7fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "921ba614",
   "metadata": {},
   "source": [
    "#### Ein Teilschritt der Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91228f13",
   "metadata": {},
   "source": [
    "Wir schreiben uns als erstes eine Funktion, welche den aktuellen Zustand der Population (beschrieben durch $S, I, R$) sowie die Parameter $\\beta$ und $\\gamma$ entgegennimmt und uns die Veränderung von $R$ für einen Tag nach obiger Formel berechnet und zurück gibt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc18c2f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nextDayR(S, I, R, beta, gamma):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4398231c",
   "metadata": {},
   "source": [
    "#### Ein Schritt der Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ff812ff",
   "metadata": {},
   "source": [
    "Nun machen wir das ganze für alle drei Grössen $S, I, R$ zusammen. Der Rückgabewert der Funktion ist nun nicht mehr nur eine Zahl, sondern ein Tupel mit den neuen Werten für $S, I$ und $R$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0731b64c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Differentialgleichungen des S.I.R. Modells\n",
    "def nextDay(S, I, R, beta, gamma):\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54f4f3f",
   "metadata": {},
   "source": [
    "Mit Hilfe der Funktion `nextDay` können wir nun den Verlauf einfach simulieren. Wir initialisieren die Parameter und rufen dann die Funktion ```nextDay``` innerhalb einer Schleife auf. In jedem Schleifendurchlauf wird ein Tag simuliert. \n",
    "Dabei speichern wir den jeweiligen Zustand der Population in Listen, damit wir diesen am Ende visualisieren können. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "89eab3f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6V0lEQVR4nO3deXxU5d338c9vJntCNhKWECAgAUwCSSAgiAiICyqitlaoivuNdam23m3Vx8fb5ZHWu1rrhlpbrRuKiKLUFUVBlDVAEMIiW4RASCAhZF/nev6YCQRISICZOTPh93695jUzZ845881R8su5rnOuS4wxKKWUUsdjszqAUkop36fFQimlVJu0WCillGqTFgullFJt0mKhlFKqTQFWB/CUuLg4k5SUZHUMpZTyK6tWrdpvjIk/enmHLRZJSUlkZ2dbHUMppfyKiPzc0nJthlJKKdUmLRZKKaXapMVCKaVUm7RYKKWUapMWC6WUUm3yWLEQkddEpEhE1jdbFisiX4nIFtdzTLPPHhCRrSKyWUQuarZ8qIisc332nIiIpzIrpZRqmSfPLF4HJhy17H5ggTEmGVjgeo+IpABTgFTXNi+KiN21zUvANCDZ9Th6n0oppTzMY/dZGGO+E5GkoxZfDox1vX4DWAjc51o+yxhTC+wQka3AcBHJAyKNMUsBRORN4Argc0/lfv2HHVTVNxIdGkR0WCDRoYFEhQXSpVMIcRFB6ImNUup05O2b8roaYwoAjDEFItLFtbwHsKzZevmuZfWu10cvb5GITMN5FkKvXr1OKuA7K3byU2FFi5+FBNpIjAmjZ0wofeMjGJwYxeDEaJI6h2kRUUp1aL5yB3dLv2nNcZa3yBjzCvAKQFZW1knN6jT/92OoqW+ktKqe0uo653NVHYVltewqqWLXgSp2lVSzdHsxr37vACA6LJDRyfGcf2YXxg3sQmRI4Ml8tVJK+SxvF4tCEenuOqvoDhS5lucDPZutlwjscS1PbGG5R4UE2ukWZadbVEir6zQ0OvipsIIf80vJ/vkA324q4j9r9xAaaOeKzASmjkgiJSHS01GVUsorvF0s5gE3AE+4nj9utvwdEXkaSMDZkb3CGNMoIuUiMgJYDlwPPO/lzC0KsNtISYgkJSGSKcN70egw5Ow6wOyV+cxds5t3V+zigpSu3DdhAP26dLI6rlJKnRLx1BzcIvIuzs7sOKAQeBj4CJgN9AJ2Ar8yxpS41n8QuBloAH5njPnctTwL55VVoTg7tn9r2hE6KyvLWDWQ4MGqet5alsfLi7ZTU9/Ib89L5o5xZxBo19talFK+TURWGWOyjlnuqWJhNSuLRZPiiloe/c8G5q3dw5Be0bw8dShdOrXetKWUUlZrrVjon7oe1DkimOd+ncnzv85kY0E5l7/wAz8VllsdSymlTpgWCy+4LD2BObePpNFh+PUry7RgKKX8jhYLL0lNiGLWtBHYbcI1/1zG7tJqqyMppVS7abHwor7xEbzzX2dRW+/gv97IpqquwepISinVLlosvKxfl048d00mG/eW8eDc9W1voJRSPkCLhQXGDejC3eclM3fNbr5YX2B1HKWUapMWC4vcdV4/0npE8uDc9RRX1FodRymljkuLhUUC7Tb+9qsMDlbX89T8n6yOo5RSx6XFwkIDunVi6sjevLdyJ5v2llkdRymlWqXFwmL3jE+mU0gg0z/daHUUpZRqlRYLi0WHBfHb8/qxeMt+svNKrI6jlFIt0mLhA645qxcxYYG8tHCb1VGUUqpFWix8QFhQADeN6sOCTUXad6GU8klaLHzEDSOTCA+y849F262OopRSx9Bi4SOiwgK5amgin64r4EBlndVxlFLqCFosfMivz+pFXYODD1bnWx1FKaWO4O1pVdVxDOwWSWavaN5dsZNbzumDiFgd6bTQ4GigvK6c6oZqahpqqG6sprq+mprGGuf7BufrusY6GhwN1DvqaXA0HPkwh1/XO+ppNI04HA4MBodxPhtz1Gsch54xtPyZcQCHP2tJa8tbX9zKflqZCO1Ev7ejTqjmT96+5G2C7EFu3acWCx/z6+G9+NOcH1mZd4DhfWKtjuO3quqrKKgsoLCqkKKqokOPkpoSymrLKKsr42DtQcrqyqiorzjp7wmwBRBoCyRAAgiwHX7YxY7dZkcQbGJDEEScDxs2bOI8qW/6zCY2EI75zG6zY+PwZ61p7Q8LoZU/OFpdfGL7OeHvbeMz5R6eOMZaLHzMxMHdefjjXD7K2a3Foh0q6irYWLKRjcUb2X5wOzvLd/LzwZ8pqi46Zt2o4ChiQ2KJCoqiS1gX+kX3Iyo4isjgSCKDIgkLCCM0IJSQgBDnwx5y6H1oQChBtqBjCoKe/anThRYLHxMWFMD5KV35Yv1eHp2USqBdu5Wayy/PZ1nBMlbuXcmG4g3kleUd+iw2JJZenXoxMmEkSVFJJIQn0CWsC13DuhIfFk9IgM5/rtTJ0mLhgyYO7s5/1u5hybZixvSPtzqOpRzGwdp9a/ky70sW7lrI7ordAMSFxjEobhAT+04kpXMKKZ1T6Bza2dqwSnVgWix80Jj+8XQKDuCTtXtO22KRdzCPOT/N4fO8zymqKiLIFsTZCWczNWUqI7qPoG9UX20CUsqLtFj4oJBAOxemduOL3L1Mv3IQQQGnR1OUMYale5byeu7rLC1YSoAEcE7iOdw79F7G9hxLeGC41RGVOm1psfBRF6d144PV+azMK2FUvzir43jckt1LmLF2Bj/u+5EuYV24K+Muftn/l8SFdvyfXSl/oMXCR43qF0dwgI0FG4s6dLHYVbaLv678KwvzF9I9vDsPjXiIK/pd4fZrxJVSp0aLhY8KDbIz8ozOLNhUyEMTz+xw7fONjkZez32dF3NeJMAWwL1D7+W6M68j0B5odTSlVAtOj8ZwPzV+YBd+Lq5i+/5Kq6O4VUFFAbfOv5VnVj/DuYnnMu+KedyUdpMWCqV8mJ5Z+LBxA7vAx7l8s7GIM+IjrI7jFisKVvD7hb+nwdHA/xv1/7j8jMs73FmTUh2Rnln4sMSYMAZ07cS3m4+9G9kfzflpDrd9dRvxofG8f9n7XNHvCi0USvkJLRY+bnRyHNk/H6CmvtHqKCfNGMOMnBk8uvRRzko4i7cueYtekb2sjqWUOgGWFAsR+b2I5IrIehF5V0RCRCRWRL4SkS2u55hm6z8gIltFZLOIXGRFZquc3a8zdQ0OVv18wOooJ6WpULy89mWu6HcFL5z3Ap2COlkdSyl1grxeLESkB3A3kGWMSQPswBTgfmCBMSYZWOB6j4ikuD5PBSYAL4qI3du5rTK8T2cCbMIPW/dbHeWkzMiZwT9+/AdX9ruSR89+lACbdpMp5Y+saoYKAEJFJAAIA/YAlwNvuD5/A7jC9fpyYJYxptYYswPYCgz3blzrRAQHkN4zmh+2FVsd5YS9/9P7/OPHf3BFvyt45OxHDg29rZTyP17/12uM2Q08BewECoCDxpj5QFdjTIFrnQKgi2uTHsCuZrvIdy07bYw6ozPr8ks5WF1vdZR2W1awjD8v+zOjeozi4ZEPa6FQys9Z0QwVg/NsoQ+QAISLyHXH26SFZS1OxSUi00QkW0Sy9+3bd+phfcTZ/eJwGFi+3T/OLvIO5nHvwntJikriyXOf1KYnpToAK/7cOx/YYYzZZ4ypBz4EzgYKRaQ7gOu56XrRfKBns+0TcTZbHcMY84oxJssYkxUf33FGa83sFU1QgI2VeSVWR2lTXWMdf1j0B+xi5/nzntfObKU6CCuKxU5ghIiEifMi+/HARmAecINrnRuAj12v5wFTRCRYRPoAycAKL2e2VHCAnfTEKFbm+f4VUc+ufpbNBzbz+KjHSeyUaHUcpZSbeL19wBizXETmAKuBBmAN8AoQAcwWkVtwFpRfudbPFZHZwAbX+ncaY/z3poOTlJUUyz+/2051XSOhQb55MdiS3Ut4c8ObTB4wmTE9x1gdRynlRmJMi83/fi8rK8tkZ2dbHcNtvtlUyM2vZ/Puf41g5Bm+NyNcWV0ZV3x0BZFBkcyaOEunMFXKT4nIKmNM1tHL9RIVPzG0VywAq372zX6LF9a8QHFNMX8e/WctFEp1QFos/ERUWCADunbyyX6LDcUbeG/ze0weMJmUzilWx1FKeYAWCz+SlRTD6p8P0OjwnaZDh3Hw+LLHiQmO4a7Mu6yOo5TyEC0WfmRo7xjKaxvYWlRhdZRD5m6Zy7r96/jDsD8QGRRpdRyllIdosfAj6T2jAVi7q9TSHE1qGmp4MedFMuIzuLTPpVbHUUp5kBYLP9KnczidQgJYm19qdRQA3tv8HkXVRdw95G6dl0KpDk6LhR+x2YTBiVE+USwq6yt5dd2rjOw+kmHdhlkdRynlYVos/Ex6YjSbCsotnwzp7Q1vc6D2AL/N/K2lOZRS3qHFws+k94ymwWHI3VNmWYbyunLeyH2DcT3HMSh+kGU5lFLeo8XCz2T4QCf3h1s+pLy+nNvSb7Msg1LKu7RY+JmukSF0jQzmR4v6LRocDczcOJOsrlmkdk61JINSyvu0WPih9MRo1uYftOS7v/75awoqC7g+5XpLvl8pZQ0tFn4ovWc0O/ZXUlpV59XvNcbwRu4bJEUm6aiySp1mtFj4oaZ+ix+9fHaxpmgN64vXMzVlqk6TqtRpRv/F+6G0HlGA9zu539n0DlHBUVx2xmVe/V6llPW0WPihqNBA+sSFs263984sDtQc4Jud33BZ38sIDQj12vcqpXyDFgs/lZIQ6dV7LT7Z/gn1jnquTL7Sa9+plPIdWiz8VGpCJLtLq73SyW2M4cMtHzIobhD9Y/p7/PuUUr5Hi4WfSk1w9lts8MLZxfr969laulXPKpQ6jWmx8FOpCc65I7zRFPXBlg8IDQjl4qSLPf5dSinfpMXCT8VFBNM1MpjcPZ7t5K5uqOaLvC+4oPcFRARFePS7lFK+S4uFH0tNiPL4mcWi/EVU1ldy+RmXe/R7lFK+TYuFH0tNiGTbvgqq6zw3XPmXO74kLjSOoV2Heuw7lFK+T4uFH0tNiMRhYNNez5xdVNZXsnj3Yi7ofQF2m90j36GU8g9aLPxY0xVRnmqK+nbXt9Q21nJxH+3YVup0p8XCjyXGhBIZEsCGAs8Uiy93fEnXsK6kx6d7ZP9KKf+hxcKPiYjH7uQ+WHuQ7/d8z0VJF+mggUopLRb+LjUhik0FZTQ0Oty63292fkODo4EJSRPcul+llH/SYuHnUhMiqW1wsH1/pVv3u2DnAhLCE0iLS3PrfpVS/kmLhZ873MntvpvzqhuqWVawjHG9xiEibtuvUsp/BVgdQJ2aM+LDCQ6wkbu7jCsz3bPPZXuWUdtYy5hEnQ1PeUd9fT35+fnU1NRYHeW0ERISQmJiIoGBge1a35JiISLRwL+ANMAANwObgfeAJCAPuNoYc8C1/gPALUAjcLcx5kuvh/ZRAXYbA7t1cmsn96L8RUQERpDVNctt+1TqePLz8+nUqRNJSUl6NusFxhiKi4vJz8+nT58+7drGqmaoZ4EvjDEDgXRgI3A/sMAYkwwscL1HRFKAKUAqMAF4UUT0DrFmUhKiyN1zEGPMKe/LYRws3LWQUT1GEWhv318cSp2qmpoaOnfurIXCS0SEzp07n9CZnNeLhYhEAucCrwIYY+qMMaXA5cAbrtXeAK5wvb4cmGWMqTXG7AC2AsO9mdnXpSZEUlbTQP6B6lPeV+7+XIprirUJSnmdFgrvOtHjbcWZRV9gH/BvEVkjIv8SkXCgqzGmAMD13MW1fg9gV7Pt813LjiEi00QkW0Sy9+3b57mfwMekuIYrd8fNeQvzF2ITG6N7jD7lfSmlOg4rikUAMAR4yRiTCVTianJqRUvlr8X2FmPMK8aYLGNMVnx8/Kkn9RNndovEJu4Z9mPRrkVkdskkOiT61IMp5UemT59OamoqgwcPJiMjg+XLl1uSIycnh88+++zQ+3nz5vHEE08AcOONNzJnzpxjtlm4cCETJ070aC4rOrjzgXxjTNN/iTk4i0WhiHQ3xhSISHegqNn6PZttnwjs8VpaPxAaZKdPXPgpz5q3t3Ivmw9s5t6h97opmVL+YenSpXzyySesXr2a4OBg9u/fT12d56csbklOTg7Z2dlccsklAEyaNIlJkyZZkqW5Ez6zEJEYERl8sl9ojNkL7BKRAa5F44ENwDzgBteyG4CPXa/nAVNEJFhE+gDJwIqT/f6OKjUhig2neK/F0j1LATinxznuiKSU3ygoKCAuLo7g4GAA4uLiSEhIICkpif379wOQnZ3N2LFjAVi0aBEZGRlkZGSQmZlJeXk5AH/9618ZNGgQ6enp3H+/s8Fk27ZtTJgwgaFDhzJ69Gg2bdoEOM8SfvOb3zB69Gj69+/PJ598Ql1dHf/zP//De++9R0ZGBu+99x6vv/46d91116GsX3/99RHbHK2yspKbb76ZYcOGkZmZyccff3zMOiejXWcWIrIQmORaPwfYJyKLjDEn+yfob4GZIhIEbAduwlm4ZovILcBO4FcAxphcEZmNs6A0AHcaYzw3gYOfSk2IZN7aPRyorCMmPOik9rFkzxLiQ+PpF93PzemUar9H/5Pr9rnlUxIiefiy1FY/v/DCC3nsscfo378/559/PpMnT2bMmNYv8njqqaeYMWMGo0aNoqKigpCQED7//HM++ugjli9fTlhYGCUlJQBMmzaNl19+meTkZJYvX84dd9zBN998A0BeXh6LFi1i27ZtjBs3jq1bt/LYY4+RnZ3NCy+8AMDrr79+xHe3tE1z06dP57zzzuO1116jtLSU4cOHc/755xMeHn4yh+6Q9jZDRRljykTkVuDfxpiHReTHk/1SY0wO0NJF/ONbWX86MP1kv+900LyTe1S/uBPe3mEcLCtYxrmJ5+pVKeq0ExERwapVq1i8eDHffvstkydPPtRP0JJRo0Zx7733cu211/KLX/yCxMREvv76a2666SbCwsIAiI2NpaKigiVLlvCrX/3q0La1tbWHXl999dXYbDaSk5Pp27fvobOO42lrm/nz5zNv3jyeeuopwHlZ8s6dOznzzDNP6Jgcrb3FIsDVj3A18OApfaPyiObDfpxMsdhYspHS2lJGJox0dzSlTsjxzgA8yW63M3bsWMaOHcugQYN44403CAgIwOFwDtLZ/J6E+++/n0svvZTPPvuMESNG8PXXX2OMOeYPLYfDQXR0NDk5OS1+59Hrt+cPtba2McbwwQcfMGDAANypvX0WjwFfAluNMStFpC+wxa1J1CmJDQ+ie1TISZ++N/VXjOg+wp2xlPILmzdvZsuWw7/ScnJy6N27N0lJSaxatQqADz744NDn27ZtY9CgQdx3331kZWWxadMmLrzwQl577TWqqqoAKCkpITIykj59+vD+++8Dzl/ka9euPbSf999/H4fDwbZt29i+fTsDBgygU6dOh/pAWtLSNs1ddNFFPP/884du0l2zZs0pHh2ndhULY8z7xpjBxpg7XO+3G2N+6ZYEym1Sup/83BZL9ixhYOxA4kJP/KxEKX9XUVHBDTfcQEpKCoMHD2bDhg088sgjPPzww9xzzz2MHj0au/3wwBHPPPMMaWlppKenExoaysUXX8yECROYNGkSWVlZZGRkHGoGmjlzJq+++irp6emkpqYe0eE8YMAAxowZw8UXX8zLL79MSEgI48aNY8OGDYc6uI/W0jbNPfTQQ9TX1zN48GDS0tJ46KGH3HKMpD1DRIhIPPBfOMdtOtR0ZYy52S0pPCArK8tkZ2dbHcOrnp6/mRe+3UruoxMIDWr/iChV9VWMmjWKqSlT9bJZZYmNGzeecpu6v7nxxhuZOHEiV111lWUZWjruIrLKGHNMn3J7+yw+BhYDX+MczE/5oJSEKBwGNheWk9Ezut3bZRdm0+BoYGR37a9QSrWsvcUizBhzn0eTqFOW6roiKnfPwRMqFkv3LCXYHsyQrkM8lEwpdbSjL4n1de3t4P5ERC7xaBJ1yhJjQokMCTjhTu4Ve1eQ0SWDYHuwh5Ippfxde4vFPTgLRo2IlLse7r1rRp0yESEl4cQ6uUtrSvnpwE8M76YD+SqlWtfeq6E6GWNsxpgQ1+tOxphIT4dTJy6lexSb9pbR6Gjf3BarCp2XBQ7rNsyTsZRSfq7dAwmKyCSc81AALDTGHDsoibJcakIkNfUOduyvoF+XTm2uv7JwJaEBoaR1TvNCOqWUv2rXmYWIPIGzKWqD63GPa5nyMSmHOrnb1xS1cu9KMuIzdFY8ddqLiIhoc53FixeTmppKRkYG1dUnNtnYRx99xIYNGzySyxva22dxCXCBMeY1Y8xrOKc31Q5vH9SvSwRBAbZ2dXI39VdoE5RS7TNz5kz+8Ic/kJOTQ2ho6Alte7LFwlecyBDl0c1eR7k5h3KTQLuNAV07tevMQvsrlDrWwoULGTt2LFdddRUDBw7k2muvxRjDv/71L2bPns1jjz3GtddeC8CTTz7JsGHDGDx4MA8//PChfbz55psMHjyY9PR0pk6dypIlS5g3bx5//OMfycjIYNu2ba0OXb5jxw5GjhzJsGHD3Hb3tTu0t8/iL8AaEfkW58x15wIPeCyVOiUp3SOZv2FviwObNbdi7wpCA0JJ7WzNwG1Ktejz+2HvOvfus9sguLj9Ledr1qwhNzeXhIQERo0axQ8//MCtt97K999/f+iu6/nz57NlyxZWrFiBMYZJkybx3Xff0blzZ6ZPn84PP/xAXFwcJSUlxMbGMmnSpCPu2B4/fnyLQ5ffc8893H777Vx//fXMmDHDvcfhFLSrWBhj3nXNaTEMZ7G4zzWJkfJBqT0ieS97F3vLauge1fqp8spC7a9QqiXDhw8nMTERgIyMDPLy8jjnnCMnBZs/fz7z588nMzMTcI4vtWXLFtauXctVV11FXJxznLXY2Nhj9n+8oct/+OGHQ4MWTp06lfvu8437ods7+dEoIMcYM09ErgP+JCLPGmN+9mw8dTIO3cm9u6zVYlFaU8qWA1u4OPNib0ZTqm0ncAbgKU0z5oFz6PKGhoZj1jHG8MADD3Dbbbcdsfy5555rc6jxEx263Be0t8/iJaBKRNKBPwI/A296LJU6JQO7RSLinAipNWuKnMMW6xAfSp2ciy66iNdee42KigoAdu/eTVFREePHj2f27NkUFxcDHJoxr/nQ48cbunzUqFHMmjULcHao+4r2FosG4xye9nLgOWPMs0DbF/ErS4QHB9Cnczi5x5mTe03RGgJtgaTF6f0VSp2MCy+8kGuuuYaRI0cyaNAgrrrqKsrLy0lNTeXBBx9kzJgxpKenc++9zpGcp0yZwpNPPklmZibbtm1rdejyZ599lhkzZjBs2DAOHmz937C3tXeI8kXAFzjnyj4X2IezWWqQZ+OdvNNxiPLm7npnNTm7Svn+vvNa/HzqZ1MxGN6+5G0vJ1PqWKfjEOW+4ESGKG/vmcVkoBa4xdWx3QN48lSDKs9J6xFF/oFqDlTWHfNZbWMtucW5DOmiTVBKqfZps1iIiB142xjztDFmMYAxZqcxRvssfNjgROetMGvzS4/5LHd/LvWOejK6ZHg3lFLKb7VZLIwxjTg7t/VGPD8yqEcUIrB217Ftnk2d21oslFLt1d6b8mqAdSLyFVDZtNAYc7dHUqlT1ikkkH7xEfzYwplFTlEOSZFJxIYce/23Ukq1pL3F4lPXQ/mR9J7RLNxcdMSd3A7jYM2+NZzXs+WOb6WUakl77+B+Q0RCgV7GmM0ezqTcJL1nNHNW5bO7tJrEmDAA8g7mcbD2IJldMi1Op5TyJ+0dovwyIAfn5bOISIaIzPNgLuUGGYnRwJH9FtpfoVTL7HY7GRkZpKWlcdlll1FaWmp1pBP2yCOP8NRTT3lk3+29dPYRYDhQCmCMyQH6eCSRcpsB3ToRFGA74oqo1UWriQmOISkyybJcSvmi0NBQcnJyWL9+PbGxsT4ziJ8xBofDYXWME7qD++jLato3b6eyTFCAjdSESHJ2lR5allOUQ0aXDJ8ce0YpXzFy5Eh2794N0OpQ4oWFhVx55ZWkp6eTnp7OkiVLAHj66adJS0sjLS2NZ555BoD77ruPF1988dD+H3nkEf72t78BLQ9znpeXx5lnnskdd9zBkCFD2LVrV6vDoU+fPp0BAwZw/vnns3mz53oJ2tvBvV5ErgHsIpIM3A0s8Vgq5TbpidG8t3IXDY0OSutK2Fm+k6v6X2V1LKVa9b8r/pdNJZvcus+BsQO5b3j7Rm9tbGxkwYIF3HLLLQBMmzatxaHE7777bsaMGcPcuXNpbGykoqKCVatW8e9//5vly5djjOGss85izJgxTJkyhd/97nfccccdAMyePZsvvvii1WHOe/XqxebNm/n3v//Niy++2Op64eHhzJo1izVr1tDQ0MCQIUMYOnSoW49dk/YWi98CD+K8i/sd4EvgcY8kUm6V3jOK15fksXVfBfm1OQDaua1UC6qrqw8NRz506FAuuOCC4w4l/s033/Dmm857k+12O1FRUXz//fdceeWVhIeHA/CLX/yCxYsXc/fdd1NUVMSePXvYt28fMTEx9OrVi+eee67FYc579epF7969GTFiBND6cOjl5eVceeWVhIU5L2CZNGmSx47PcYuFiIQAvwH6AeuAkcaYY8fqVT4r3dXJnbOzlF2yhiBbECmdU6wNpdRxtPcMwN2a+iwOHjzIxIkTmTFjBjfeeONxhxI/2vHG2rvqqquYM2cOe/fuZcqUKYfWb2mY87y8vEMF53jrPfPMM15rUm6rz+INIAtnobgYcFs3u4jYRWSNiHzieh8rIl+JyBbXc0yzdR8Qka0isllELnJXhtNBUudwosMCWbOzlJyiHNLi0giyB1kdSymfFRUVxXPPPcdTTz1FaGhoq0OJjx8/npdeeglwNl2VlZVx7rnn8tFHH1FVVUVlZSVz585l9OjRgHPU2VmzZjFnzpxDs+W1Nsz50Vpb79xzz2Xu3LlUV1dTXl7Of/7zH48dl7aaoVKaRpYVkVeBFW787nuAjUCk6/39wAJjzBMicr/r/X0ikgJMAVKBBOBrEenvGoZEtcFmE4b2imHFz3spjd/A9anXWx1JKZ+XmZlJeno6s2bNYubMmdx+++08/vjj1NfXM2XKFNLT03n22WeZNm0ar776Kna7nZdeeomRI0dy4403Mnz4cABuvfXWQ01HqamplJeX06NHD7p37w44hznfuHEjI0eOBCAiIoK3334bu91+RJ7W1hsyZAiTJ08mIyOD3r17HypMnnDcIcpFZLUxZkhr70/6S0UScZ61TAfuNcZMFJHNwFhjTIGIdAcWGmMGiMgDAMaYv7i2/RJ4xBiz9HjfcboPUd7cSwu38dR3nxHW+xWeP+95xvYca3UkpY6gQ5Rbw51DlKeLSJnrUQ4MbnotIq1Pw9a2Z4A/Ac0vHu5qjCkAcD13cS3vAexqtl6+a9kxRGSaiGSLSPa+fftOIV7HMiwpBnuocwbcjPgMa8MopfzScYuFMcZujIl0PToZYwKavY483ratEZGJQJExZlV7N2kpWit5XzHGZBljsuLj408mXoc0KDGKgPA8Iu09iA6JtjqOUsoPtffSWXcaBUwSkUuAECBSRN4GCkWke7NmqKZennygZ7PtE4E9Xk3s5wLtQmDYLqjRS2aV72o+4KXyvPbMktpce+/gdhtjzAPGmERjTBLOjutvjDHXAfOAG1yr3QB87Ho9D5giIsEi0gdIxr0d7R3ettJtOKSK4uIEquv0ugDle0JCQiguLj7hX2Dq5BhjKC4uJiQkpN3bWHFm0ZongNkicguwE/gVgDEmV0RmAxuABuBOvRLqxDQNHlhb0ZvVOw8wql+cxYmUOlJiYiL5+floX6P3hISEkJiY2O71LS0WxpiFwELX62JgfCvrTcd55ZQ6CTlFOUQHx1DVGMfSbcVaLJTPCQwMpE8fHZvUl3m9GUp535qiNQzpksngxGiWbNtvdRyllB/SYtHB7a/eT35FPpldMjn7jM6szT9IRa2O2KKUOjFaLDq4nKIcwDnZ0dlnxNHoMKzcUWJtKKWU39Fi0cGtKTo8eODQ3jEE2W3aFKWUOmG+dDWU8oCcohxS41KdgwfaYUjvaH7YWmx1LKWUn9Eziw6spqGGDSUbjphve3RyPBsKyigqq7EumFLK72ix6MDW719Pg6OBzPjDd26PG+AccmvhT3o9u1Kq/bRYdGA5+3IAjjizOLN7J7pGBrNw87Fj5iulVGu0WHRgOUU5JEUmERNyaB4pRIRxA7qw+Kf91Dc6jrO1UkodpsWig3IYBzn7co44q2gydkAXymsbWPXzAe8HU0r5JS0WHVTewTwO1h4ks8uxI82ekxxHoF1YsLHQgmRKKX+kxaKDaho8sKUzi4jgAEb1i+OL3L06yqdSql20WHRQOftyiA6Opk9ky4OzXZLWnV0l1eTuOZUJD5VSpwstFh1UTlEOGfEZrU4mc0FKV+w24fP1BV5OppTyR1osOqCSmhLyyvJabIJqEhMexMi+nflsnTZFKaXapsWiA2oaPLClzu3mLh7UjR37K9lYUO6FVEopf6bFogNaVbiKIFsQqXGpx13vkrTuBNqFD1fneymZUspfabHogLILsxkUP4hge/Bx14sJD2L8wK58lLNbb9BTSh2XFosOpryunE0lm8jqmtWu9X85NJH9FXV8p2NFKaWOQ4tFB7OmaA0O4yCrW/uKxdgB8XQOD+L9bG2KUkq1TotFB5NdmE2ALYD0+PR2rR9ot3FlZg++3lhIoQ5brpRqhRaLDmbV3lWkdU4jNCC03dtcN6I3jcYwc/lODyZTSvkzLRYdSFV9FbnFue1ugmqSFBfOuAFdeGf5TuoatKNbKXUsLRYdSE5RDo2msd2d283dcHYS+ytq+XTdHg8kU0r5Oy0WHUh2YTZ2sR/3zu3WjO4XR3KXCF5euB2HQ+/oVkodSYtFB5JdmE1K5xTCA8NPeFubTbjrvH5sLizny9y9HkinlPJnWiw6iMr6StbtX8ewbsNOeh8TByfQNy6cZxds0bMLpdQRtFh0ENl7s2lwNDAyYeSRH9QchJ+XQN73UH38mfHsNuHOcf3YtLecz3Q0WqVUMwFWB1DusbRgKSH2kMODB1aVwILHIOcdaKx1LhM7pEyC8Q9DbMvzXFyR2YN/Lt7OE59v4vwzuxISaPfST6CU8mV6ZtFBLNmzhKFdhzrHg9q/BV4ZA2vegsxr4doP4LoPYcTt8NN8eHk0bP6ixf3YbcL/TEwh/0A1r/2ww8s/hVLKV3m9WIhITxH5VkQ2ikiuiNzjWh4rIl+JyBbXc0yzbR4Qka0isllELvJ2Zl+3t3IvOw7ucDZBle2Bt66E+mq4eT5M/Dsknw/9xsNF0+HO5dC5L7w7BVa+2uL+zu4XxwUpXXnhm63sKqny8k+jlPJFVpxZNAD/bYw5ExgB3CkiKcD9wAJjTDKwwPUe12dTgFRgAvCiiGjbSDNL9ywFYGR8Jrx9FVSXwnUfQOLQY1eO7gk3fQH9L4JP74UfZ7e4z0cmpSLAAx+u08mRlFLeLxbGmAJjzGrX63JgI9ADuBx4w7XaG8AVrteXA7OMMbXGmB3AVmC4V0P7uKV7lhIXGkfyyregKBeufh26H2dsqKAw+NUbkDQa5v4Gti86ZpUe0aE8cMmZfL91P++s0GFAlDrdWdpnISJJQCawHOhqjCkAZ0EBurhW6wHsarZZvmtZS/ubJiLZIpK9b9/pMeS2wzhYVrCMkZ36ICtfgbN+A/3Ob3vDwBCY8g7EJcOcm6B01zGrXDO8F6OT43jsPxvYWFDmgfRKKX9hWbEQkQjgA+B3xpjj/SaSFpa12C5ijHnFGJNljMmKj493R0yft6lkEwdqDzAybxXEnuG80qm9QiJh8kxorIfZ10ND3REf22zC3ydnEBUayB0zV1NWU+/m9Eopf2FJsRCRQJyFYqYx5kPX4kIR6e76vDtQ5FqeD/RstnkioAMYuXy/+3sARu7fCZf81dnEdCLi+sHlM2DPalj4l2M/jgjmhWuGsKukitvfXqUDDSp1mrLiaigBXgU2GmOebvbRPOAG1+sbgI+bLZ8iIsEi0gdIBlZ4K6+vW/Tz1wyqayCu7/nta35qScokGHI9fP935817RxneJ5b//eVgfthazJ/mrKVR7+5W6rRjxZnFKGAqcJ6I5LgelwBPABeIyBbgAtd7jDG5wGxgA/AFcKcxptGC3D5nf/V+fizZyNjKSrjw8VPb2UV/cd6o9+FtLd7p/cuhifzxogF8lLOHP76vBUOp043X7+A2xnxPy/0QAONb2WY6MN1jofzUoi3zABjTcxzEDzi1nQVHwC//Ba9eCJ/9CX75z2NWuXNcPxodhqe/+onq+kaevjqD0CC9ilmp04Hewe3HFua+Q0JDA/3HPOSeHfYYCuf+CdbNhty5La5y9/hk/u+lZ/JF7l6mvLKUPaXV7vlupZRP02Lhp6oO5LG0Zi9jw3ohcf3ct+PR90LCEPjk91De8lDlt47uyz+uG8rWogoufnYxn/6ogw4q1dFpsfBTixc9Qq1NGD/kdvfu2B4Iv3jFOVzIx3dBK3dvX5jajU/vHk1SXDh3vrOa381aQ1FZjXuzKKV8hhYLf1ReyJcFPxAngQxNvsz9+49Lhgseg61fwarXW10tKS6cOb8Zyd3n9eOzdXsZ99RCZny7laq6BvdnUkpZSouFH6pc/De+Cwnmgt4XYLd5qIN52H9BnzHw5YNQsr3V1QLtNu69cABf3Xsuo/rF8eSXmxn1xDc8+/UWDlTWtbqdUsq/aLHwNxVFLNr0HrU24aKBV3vue2w2uOJFsAXA3NvBcfyrlXt3DueV67OY85uRDO0dw9+//omz/rKAu95ZzaKf9tHQqDfzKeXPdPIjf7PkeT4PDaRLSOfDEx15SlQiXPIkzJ0GS56Dc37f5iZZSbH8KymWzXvLeXfFTj7K2c0nPxYQFRrI2AHxnDewC+cmxxMTHuTZ7Eopt9Ji4U8qi9m/6jW+T+jMdWdchk28cGI4+GrY/Cl8M915h3i3Qe3abEC3TjwyKZUHLhnIt5v28fXGQr7dVMTHOc6RWs6ID2do7xiG9IphQLdO9O/aifBg/d9RKV+l/zr9ybIZfBpsowHDlf2u9M53isClf4eflzrv7p72LQQEt3vz4AA7E9K6MSGtG40Ow9r8UpZtL2ZV3gHmbyhkdnb+oXV7xoaS3KUTPWNC6RETSo/oMBKiQ+gRHUpseBABdm01VcoqWiz8RVUJZvkrzE3sTnrcAPpG9/Xed4d3hstfgHeuhm+nO6+UOgl2mzCkl/NsAsDhMOwsqeKnwnJ+Kixnc2EFWwrLWbmjhPLaY6+oigoNJDY8iNjwIGLCgogJCyQ8OICI4ADCgu2EBwUQHhxAeJCdsOAAQgPtBNqFQLuNoAAbQXYbgQE2Au1CkGtZoN1GgE1wDlmmlGqNFgt/sfxlfpRatpkaHvHWWUVz/S+CoTfCD885J01KvuCUd2mzCUlx4STFhXNharcjPiurqWf3gWr2lDofxZV1HKiscz5X1bG7tJrcPQepqG2gsraBUx2qKsAm2EQQcRa1ptc2Edd7EHE+28VZXGw25+umdZsXnKZXzWuQuJa2VJfa3FaO3MeRyzhmRTl2kTqNzJo2guAA914pqcXCH1QfgGUvMyvxTMJstUzoM8GaHBf9BfJXwQe3wm2LICbJY18VGRJIZPdAzuwe2ea6xhhqGxxU1jZQVddIZZ2zgNTWO6htdFDf4KC+0VDX2Eh9g6Gu0UFdg4P6RuejrsFBozE4jPNsx9H02hjXe9drY3A4OPLzZtsdyuOabqX5/YxNr02zqVgOLzt2Pdpcz7S6bUvrqdOLtDr83snTYuEPFj9NUUMFX5gyJidPITww3JocQWEw+U34x1jnZEk3z3fOuGcxESEk0E5IoJ3OVodRqoPSHkNfV7oLlv+DWf3OotE4uHbgtdbmie0Lv/gHFKyFea0PB6KU6li0WPi6b/9MtcD7lDG251h6RvZsextPG3AxnPcQrHsfvv2z1WmUUl6gxcKX7V0Ha99lVsp5lNaVcVPaTVYnOmz0f0PmdfDdX2HNTKvTKKU8TPssfJUx8OX/oTI0itdq8zk74WzP37F9IkRg4jPOZrL/3A2h0TDwUqtTKaU8RM8sfNW692HHd8wcNIHSuoPclXGX1YmOZQ+EyW9D9wyYfQP8NN/qREopD9Fi4YuqD8CX/4eiHhm8eiCHsT3HMii+fcNseF1IJFz3AXRNgfeugy1fWZ1IKeUBWix80Wd/hKoSnuo1gAZHA3/K+pPViY4vNBqmfuScB/zdKbB2ltWJlFJupsXC16ybA+veZ9nw6/l871JuHnSzb1wB1ZawWLjxU+g9CubeBouf1stqlepAtFj4kuJt8Om9HEwcyv+tWE/vyN7cknaL1anaLyQSrp0DaVfBgkdhzs1QW2F1KqWUG2ix8BW15TDrWozYeLRXf4qri3li9BOEBFh/h/QJCQiCX/4Lzn8ENnwE/zwPijZanUopdYq0WPiCxnqYcwvs38xrI67lq4IfuDPzTtLi0qxOdnJEnBMlTZ0LVcXwj3OdzVKNOje3Uv5Ki4XVHA746HbY8iVfnjONZ/LmMSFpAjen3Wx1slPXdyzcsQz6T3A2S716PuxaaXUqpdRJ0GJhpYZa+OBmWPc+X4+4kfv3zCcjPoPHz3ncO7PgeUNEPFz9Jlz1bygrcBaMD26FAz9bnUwpdQI6yG8kP1SxD97+JeTO5cOzruUPRYtI65zGi+e/SLC9/TPR+QURSPsF/HYVjP4DbJgHzw+Bj+6E/VutTqeUagfpqOPdZ2VlmezsbKtjtGz7IvhwGjU1B3gy42JmF69mZPeRPDPuGcICw6xO53kHd8OS52DV69BQA33HOSdWGnCJs4NcKWUZEVlljMk6ZrkWCy8q3wtfPwpr32Fll748Gh/Hz1V7uSn1Ju4ecjcBttNsqK6KIsh+DVa/BWX5EB4PKZfDmZc579ewB1qdUKnTjhYLK5Vsh5WvYla+yo8B8ErSIL6r3UuPiB48evajnNX9LKsTWsvRCFsXwJq3nMOFNFRDSDT0G++cwjVpNHQ+Q+cHVcoLWisWp9mfsl5UsQ+2fgU/zmbfz4v4Njyc93v3ZpOjik6mkt8N+R3XnHkNoQGhVie1ns0O/S90PuqqYNs3sOkT2PYtrP/AuU5EV+ie7nx0GwzdB0NUT+e2SimP85tiISITgGcBO/AvY8wTFkc6rKEWSnZA4XrMrpUU5C9lfekW1gcHsTKiE+t7JQLQPyqRhwZM5tK+l1o3NaqvCwqDMyc6H8Y472rPWww7l8HeH51nIKbRua49CKJ7O2fvi+0DMX0gsruzsER0cT4H6XFWyh38ohlKROzAT8AFQD6wEvi1MWZDa9ucdDNUXRXUVUJdBdQ7X9fXHKS2tpSa6lJqKgs5WFFASVURpdX7OVC1n+L6Mnbb7ewKDGBXYCDlNudFZoESQEpcKmMSxzCm5xiSo5MRbUo5NfXVULgBCtc5m/dKtkNJnvO5vvLY9YMinH0hIVHO4UiCI52vgyOd74MiIDAUAkKcj8AQCAiFgGDX8mDncluA8yzGFtDsYQexH/Ve//sq/+bvzVDDga3GmO0AIjILuBxotVicrNvfHMl26qixCbXifDQc7xdAGAQQRUJwDD0jEkmL7U9y7AAGxQ0iOSaZILte3eNWgaGQONT5aM4YqNwPFXuhotDZeX7ouQhqy6DmIFRsc70ug7py9+cT++GiIjZAXAVEQDjq/Yk+t7K9W3K7s8h14EzgH38Q3Pad8w8dN/KXYtED2NXsfT5wTK+wiEwDpgH06tXrpL6oX/csYuvLCA4IJTgglJDAMIIDwwkJDCc4KIKQ0DhiIroRHRJDTEgM0cHRRAZF6hmD1UScNwBGxAPtnPvD0eg8i2yodXaq19c4L+Vtehx6XwuOBufDNLpeNx5e5mg86n3TugYwp/BMO9ZzBze2LnToTODWXB7l/t9H/lIsWvrJj/mvZox5BXgFnM1QJ/NF/33pqyezmfJHNruzKUop1SZ/uYM7H2g+qUMisMeiLEopddrxl2KxEkgWkT4iEgRMAeZZnEkppU4bftEMZYxpEJG7gC9xXjr7mjEm1+JYSil12vCLYgFgjPkM+MzqHEopdTryl2YopZRSFtJioZRSqk1aLJRSSrVJi4VSSqk2+cXYUCdDRPYBJzt3Zxyw341xPEVzupfmdC9/yOkPGcG7OXsbY+KPXthhi8WpEJHslgbS8jWa0700p3v5Q05/yAi+kVOboZRSSrVJi4VSSqk2abFo2StWB2gnzelemtO9/CGnP2QEH8ipfRZKKaXapGcWSiml2qTFQimlVJu0WDQjIhNEZLOIbBWR+63O05yI5InIOhHJEZFs17JYEflKRLa4nmMsyPWaiBSJyPpmy1rNJSIPuI7vZhG5yOKcj4jIbtcxzRGRS3wgZ08R+VZENopIrojc41ruU8f0ODl96piKSIiIrBCRta6cj7qW+9rxbC2n7xxPY4w+nP02dmAb0BcIAtYCKVbnapYvD4g7atlfgftdr+8H/teCXOcCQ4D1beUCUlzHNRjo4zredgtzPgL8oYV1rczZHRjiet0J+MmVx6eO6XFy+tQxxTnLZoTrdSCwHBjhg8eztZw+czz1zOKw4cBWY8x2Y0wdMAu43OJMbbkceMP1+g3gCm8HMMZ8B5Qctbi1XJcDs4wxtcaYHcBWnMfdqpytsTJngTFmtet1ObAR5xz0PnVMj5OzNVblNMaYCtfbQNfD4HvHs7WcrfF6Ti0Wh/UAdjV7n8/x/+f3NgPMF5FVIjLNtayrMaYAnP94gS6WpTtSa7l88RjfJSI/upqpmpoifCKniCQBmTj/yvTZY3pUTvCxYyoidhHJAYqAr4wxPnk8W8kJPnI8tVgcJi0s86XrikcZY4YAFwN3isi5Vgc6Cb52jF8CzgAygALgb67llucUkQjgA+B3xpiy463awjKvZW0hp88dU2NMozEmA0gEhotI2nFW97WcPnM8tVgclg/0bPY+EdhjUZZjGGP2uJ6LgLk4TzkLRaQ7gOu5yLqER2gtl08dY2NMoesfqAP4J4dP4y3NKSKBOH8BzzTGfOha7HPHtKWcvnpMXdlKgYXABHzweDZpntOXjqcWi8NWAski0kdEgoApwDyLMwEgIuEi0qnpNXAhsB5nvhtcq90AfGxNwmO0lmseMEVEgkWkD5AMrLAgH3Dol0STK3EeU7Awp4gI8Cqw0RjzdLOPfOqYtpbT146piMSLSLTrdShwPrAJ3zueLeb0qePp6V5+f3oAl+C8qmMb8KDVeZrl6ovzyoe1QG5TNqAzsADY4nqOtSDbuzhPj+tx/rVzy/FyAQ+6ju9m4GKLc74FrAN+xPmPr7sP5DwHZ3PCj0CO63GJrx3T4+T0qWMKDAbWuPKsB/7HtdzXjmdrOX3meOpwH0oppdqkzVBKKaXapMVCKaVUm7RYKKWUapMWC6WUUm3SYqGUUqpNAVYHUKojEJFGnJc4BgINOMcbesY4b6ZSyu9psVDKPaqNc6gGRKQL8A4QBTxsZSil3EWboZRyM+MckmUazgHgRESSRGSxiKx2Pc4GEJG3ROTQyMYiMlNEJolIqmtugxzXAHLJVv0sSjXRm/KUcgMRqTDGRBy17AAwECgHHMaYGtcv/neNMVkiMgb4vTHmChGJwnkXdDLwd2CZMWama+gZuzGm2qs/kFJH0WYopTynaWTQQOAFEckAGoH+AMaYRSIyw9Vs9QvgA2NMg4gsBR4UkUTgQ2PMFguyK3UEbYZSygNEpC/OwlAE/B4oBNKBLJwzMTZ5C7gWuAn4N4Ax5h1gElANfCki53kvuVIt02KhlJuJSDzwMvCCcbbzRgEFriujpuKcwrfJ68DvAIwxua7t+wLbjTHP4Rw8brDXwivVCm2GUso9Ql2znDVdOvsW0DR094vAByLyK+BboLJpI2NMoYhsBD5qtq/JwHUiUg/sBR7zeHql2qAd3EpZSETCcN6fMcQYc9DqPEq1RpuhlLKIiDRNxPO8Fgrl6/TMQimlVJv0zEIppVSbtFgopZRqkxYLpZRSbdJioZRSqk1aLJRSSrXp/wNTY6jh+HcCGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Anfangszustand\n",
    "N = 1000       # Anzahl Personen\n",
    "I = 1          # Infizierte\n",
    "R = 0          # Genesene\n",
    "S = N - I - R  # Gefährdete Personen\n",
    "\n",
    "# Parameter\n",
    "beta = 0.2        # Ansteckungsrate\n",
    "gamma = 1.0 / 12  # Genesungsrate\n",
    "\n",
    "# In diesen Listen werden die Zahlen pro Tag gespeichert\n",
    "Ss = []\n",
    "Is = []\n",
    "Rs = []\n",
    "\n",
    "# Simulation mehrerer Tage\n",
    "days = 365\n",
    "for i in range(1, days):\n",
    "    pass\n",
    "\n",
    "\n",
    "# Visualisierung\n",
    "plot(Ss, label= \"Susceptible\")\n",
    "plot(Is, label= \"Infected\")\n",
    "plot(Rs, label = \"Recovered\")\n",
    "legend()\n",
    "xlabel(\"Days\")\n",
    "ylabel(\"Persons\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b545cfcc",
   "metadata": {},
   "source": [
    "### Übungen zum selber probieren\n",
    "\n",
    "* Erweitern Sie die Simulation, so dass die Immunität wieder verloren gehen kann. Führen Sie dazu einen neuen Modellparameter $\\delta$ ein\n",
    "* Erweitern Sie die Simulation, so dass infizierte Personen zu einem bestimmten Anteil sterben können.\n",
    "* Erweitern Sie die Simulation, so dass auch neue Personen geboren werden können. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee872d36",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7c2bd0dd",
   "metadata": {},
   "source": [
    "### Referenzen\n",
    "1. [Bhaskar Krishnamachari - Simulating Covid-19 with 6 lines of code — the SIR epidemic model](https://medium.com/@bhaskark2/simulating-covid-19-in-6-lines-of-code-the-sir-epidemic-model-89aab118cb7b)\n",
    "2. [Christian Hill - Learning Scientific Programming with Python](https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}