Last active
May 26, 2020 06:26
-
-
Save stijnvanhoey/0b6a89c1ab12e4c9a4bff123cdf281b4 to your computer and use it in GitHub Desktop.
SIR model implementation tryout
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "from scipy.integrate import solve_ivp\n", | |
| "import pandas as pd\n", | |
| "import matplotlib.pyplot as plt" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## As simple and scipy-close as possible" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "# User specifies...\n", | |
| "\n", | |
| "# ... state transitions\n", | |
| "def integrate(t, y, parameters):\n", | |
| " \"\"\"Basic SIR model\"\"\"\n", | |
| " \n", | |
| " # unpacking need to be done explicitly (how to control if correct?)\n", | |
| " S, I, R = y \n", | |
| " beta, gamma = parameters\n", | |
| " \n", | |
| " # Model equations\n", | |
| " N = S + I + R\n", | |
| " dS = -beta*S*I/N\n", | |
| " dI = beta*S*I/N - gamma*I\n", | |
| " dR = gamma*I\n", | |
| "\n", | |
| " return dS, dI, dR\n", | |
| "\n", | |
| "# ... time, parameters and initial conditions\n", | |
| "time = [0, 250]\n", | |
| "parameters = {\"beta\": 0.5, \"gamma\": 0.3} # same order as definition\n", | |
| "initial_states = {\"S\": 7900000, \"I\": 10, \"R\": 0} # same order as definition\n", | |
| "\n", | |
| "# -> runs model with scipy directly.\n", | |
| "output = solve_ivp(integrate, time, list(initial_states.values()), args=[list(parameters.values())])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "[<matplotlib.lines.Line2D at 0x7f4d7674a2d0>]" | |
| ] | |
| }, | |
| "execution_count": 11, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3/8ddnZrISQoCELQn7vgsRpShVwRaXqtWKIrhV66+9va1ttb229vbetrd3ab3WeutaXMsiKtR9b13ABQmbLGFfA0jCkhDIOjOf3x9nEgIEGGAm50zyeT4c58w5Zw7vJJN3Tk7OfI+oKsYYY7zL53YAY4wxJ2ZFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHhe3ohaRJ0WkRERWRrn+ZBFZLSKrRGRWvHIZY0yikXidRy0i44GDwLOqOvQk6/YDngcuUtX9ItJJVUviEswYYxJM3PaoVfUjYF/jeSLSR0TeEpHFIjJfRAZGFn0HeEhV90eeayVtjDERzX2M+nHgB6o6GrgbeDgyvz/QX0Q+FpHPRGRSM+cyxhjPCjTXPyQiGcBXgBdEpH52SqMc/YALgDxgvogMVdWy5spnjDFe1WxFjbP3XqaqI5tYVgx8pqp1wGYRWYtT3IuaMZ8xxnhSsx36UNUDOCV8LYA4RkQWvwRcGJmfjXMoZFNzZTPGGC+L5+l5s4FPgQEiUiwitwFTgdtEZDmwCrgysvrbwF4RWQ28D/xUVffGK5sxxiSSuJ2eZ4wxJjbsnYnGGONxcfljYnZ2tvbs2TMemzbGmBZp8eLFe1Q1p6llcSnqnj17UlhYGI9NG2NMiyQiW4+3zA59GGOMx1lRG2OMx1lRG2OMx0VV1CLy48jwoytFZLaIpMY7mDHGGMdJi1pEcoEfAgWR4Ur9wPXxDmaMMcYR7aGPAJAmIgEgHdgZv0jGGGMaO2lRq+oO4D5gG7ALKFfVd45eT0TuEJFCESksLS2NfVJjjGmlTnoetYi0xxmToxdQhjNM6TRVndF4PVV9HGe8aQoKCk7rfekP/n09wVD4ZIFOZ9OHnx7F5gRBxFm3fnn90KzOfDk8v9E8n08I+AR/5L7xY3/DtA+/D/w+n7OOCEl+ISM1QGZqEm1TA2SkBJAz/DiNMS1HNG94mQhsVtVSABGZhzOu9IwTPus0PPrhRqrqQsdd3lqGJfEJtE1NIjPNKe/MxtNpxz4e2KUt+R3S3Y5tjImTaIp6G3CuiKQDVcAEIC5vO1z9m+a9sEtTA1KpgkaW6RHztOEHRcN9ZF79+uEwBMNhQqqEwkowpIRVCYadx/W34BHTYYIh5WBNkIrqOg5UBTlQXceBqjoOVAcj93Vs2VPZMP9Q7bE/zHp2TOe8ftmc3y+HsX06kpmaFJfPmTGm+Z20qFV1oYi8CCwBgsBSIoc4El1ThxcOz/LuoYdgKExFtVPo+yvrWLptP/PX72Hekh3M+Gwbfp8wIq8d5/fL4fx+2YzIzyLJb6fMG5Oo4jLMaUFBgdpYH82vNhhmybb9LFi/h/kb9rCiuIywQkZKgLF9OnJ+v2zO65tNr+w2dgzcGI8RkcWqWtDkMivqlqusspZPNu5l/vo9zF9fSvH+KgBys9Kc0o4Ud1Z6sstJjTFW1AZVZeveSuZv2MP8daV8unEvFTVB2qYE+Ovt5zAyP8vtiMa0albU5hjBUJhl28v4yfPLKa+q47k7zmVQ10y3YxnTap2oqO0vTK1UwO+joGcHZt5+DmlJfm584nM2lR50O5YxpglW1K1cfod0Ztx+DqrKtOkLKd5f6XYkY8xRrKgNfTtl8OxtY6ioCTJt+kJKKqrdjmSMacSK2gAwpFs7nr51DCUVNdw4/XP2H6p1O5IxJsKK2jQY3aM9028qYPPeQ9z81OdUVNe5HckYgxW1OcpX+mbz8A2jWL3zALc9U0hVE29XN8Y0Lytqc4yJgztz/3UjWbRlH9+dsZiaoJW1MW6yojZNumJEN/776mF8uK6UO2cvO/nws8aYuLGiNsd13dnd+dfLB/PWqi/52dwvCIdbyTizxnhMNMOcmlbstvN6cagmyP3vrqNNcoDfXDnEBnQypplZUZuT+sFFfTlUE+SxjzbRJiXAv0waYGVtTDOyojYnJSLcc8lADtYEefTDjbRNDfD9C/u6HcuYVsOK2kRFRPjtlUOprA3xh7fXkp7s59ZxvdyOZUyrEM3FbQcAcxrN6g38SlUfiFsq40k+n/CHbw2nsjbIr19dTZuUAJML8t2OZUyLd9KzPlR1raqOVNWRwGigEvhb3JMZTwr4fTw45SzO75fNPXO/4PUvdrkdyZgW71RPz5sAbFTVrfEIYxJDSsDP4zcWMLpHe+58bikbbXhUY+LqVIv6emB2UwtE5A4RKRSRwtLS0jNPZjwtLdnPQ1NHAfD8ou0upzGmZYu6qEUkGbgCeKGp5ar6uKoWqGpBTk5OrPIZD+vUNpWLBnZi7pId1Nk7F42Jm1PZo74EWKKqu+MVxiSeyQX57DlYw/trStyOYkyLdSpFPYXjHPYwrdcFA3LIaZvC84XFbkcxpsWKqqhFJB24GJgX3zgm0QT8Pq4elcv7a0vsyjDGxElURa2qlaraUVXL4x3IJJ7JBfmEwsrfluxwO4oxLZKNnmfOWJ+cDAp6tGdO4XZUbYQ9Y2LNitrExOSCfDaVHmLJtv1uRzGmxbGiNjFx6fCupCf7eX6R/VHRmFizojYxkZES4PLhXXnti50cqgm6HceYFsWK2sTM5IJ8DtWGeH2Fjf9hTCxZUZuYGd2jPb2z2/BCob2l3JhYsqI2MSMiXFuQz6It+9lkAzUZEzNW1CamrhmVi98nvLDY/qhoTKxYUZuY6pSZyoUDcpi7uJigDdRkTExYUZuYu7Ygn5KKGj5cZ8PdGhMLVtQm5i4a2InsjGSetz8qGhMTVtQm5pL8Pq4elcffi0rYc7DG7TjGJDwrahMX147OI2gDNRkTE1bUJi76dW7LWd2zeN4GajLmjFlRm7iZXJDP+pKDLNte5nYUYxKaFbWJm8uHdyUtyW9XfzHmDEV7hZcsEXlRRNaISJGIjI13MJP42qYmcemwrry6fCeVtTZQkzGnK9o96j8Bb6nqQGAEUBS/SKYlmVyQx8GaIG+u+NLtKMYkrJMWtYhkAuOBJwBUtVZV7aCjicqYXh3o2THdzqk25gxEs0fdGygFnhKRpSIyXUTaHL2SiNwhIoUiUlhaau9IM476gZoWbt7Hlj2H3I5jTEKKpqgDwCjgEVU9CzgE3HP0Sqr6uKoWqGpBTk5OjGOaRHbNqDx8Ai8str1qY05HNEVdDBSr6sLI4xdxituYqHRpl8pX++fw4uJiQmE7p9qYU3XSolbVL4HtIjIgMmsCsDquqUyLM7kgn90HavhovR0WM+ZURXvWxw+AmSLyBTAS+M/4RTIt0YRBnenQJtmu/mLMaQhEs5KqLgMK4pzFtGDJAR/fPCuXZz/dwt6DNXTMSHE7kjEJw96ZaJrN5IJ86kLKS8t2uh3FmIRiRW2azYAubRmR144XbKAmY06JFbVpVpPPzmfNlxWs2FHudhRjEoYVtWlW3xjRjZSAjzmL7I+KxkTLito0q8zIQE2vLNtJVW3I7TjGJAQratPsri3Io6ImyNurbKAmY6JhRW2a3bm9OtK9gw3UZEy0rKhNs/P5hGtH5/HJxr1s21vpdhxjPM+K2rjimtF5iMCLNlCTMScV1TsTjYm1bllpnN/PGajpzon98fvE7UjmKKpKSEOENdxwA5xpwqgqqnp4GiWsh6cb7htNOxumYbp+WcN0o3WOyHL0jEi+Y+Y1sV5UH+tpPu9oAQnQO6t3TLZ1xHZjvkVjonRdQT7fn7WEjzfsYXx/Gxr3ZA7WHuRA7QGqglVU1lVSGaw8fN9ouqqu6shldZVUBauoC9cR0hDBcJBgONgwHQqHCGrT88yp6ZjakQ+u+yDm27WiNq6ZOLgTWelJzCncbkUNVAWr2HlwJzsO7nBuFTsaposPFlNRWxHVdlL8KaQH0klPSictkNZwn+HLIOALEJAAfp+fgC+AX/wk+ZLwi79hXkACzjKf35kvfnziQ0Tw4dwL0jCvYRpx1mk0LTi/KdVPi0Qec/g3qCPWO846kRnHOGad48yLRv2/eyaSfclnvI2mWFEb16QE/Fw1MpdZC7ex/1At7dvE50XuFapK8cFiiiuKGwq5+GBxQynvrd57xPop/hS6ZXQjNyOX4TnDyc3IJSsli7SkNKeII2Xc+D41kErAZ9/WLY19RY2rJhfk8/QnW3h52Q5uGdfL7ThxsalsE69vfp03N7/J9orDfzz1i58ubbqQl5HHV/O/Sm5G7hG37LTsmOzlmcRnRW1cNbhbJsNy2/F8YXGLKupdB3fx5pY3eWPTG6zdvxaf+BjTZQw3D76Z3lm9yc3IpVN6J9v7NVGxV4lx3eSCPP715VWs3FHO0Nx2bsc5bfur9/POlnd4Y/MbLClZAsDw7OHcM+Yevt7z62SnZbuc0CSqqIpaRLYAFUAICKqqXUTAxMwVI3L57etFPF+4PeGK+lDdIf6x7R+8uflNPt35KUEN0rtdb35w1g+4pOcl5Gfmux3RtACnskd9oaruiVsS02q1S09i0pAuvLR0B7+4dBCpSX63I51QbaiWBTsW8ObmN/lg+wdUh6rp2qYrNw25iUt7XUr/9v3t2LKJKTv0YTzhurPzeWX5Tt5ZvZsrRnRzO06TVJWHlj3ErDWzqKitoH1Ke67seyWX9rqUkZ1G4hN7o6+Jj2iLWoF3RESBx1T18aNXEJE7gDsAunfvHruEplUY27sjuVlpPL9ouyeLWlX53cLfMWftHCZ2n8jV/a7m3G7nkuRLcjuaaQWi3QUYp6qjgEuA74vI+KNXUNXHVbVAVQtycuzNC+bU+HzCtQV5fLxxD9v3eWugJlXlPxf+J3PWzuHWIbdy/wX3c37e+VbSptlEVdSqujNyXwL8DRgTz1CmdfrW6DwA5i4pdjnJYarKf33+Xzy39jluGXILPx79Yzv+bJrdSYtaRNqISNv6aeBrwMp4BzOtT177dM7rm80LhcWEwu5f/FZV+Z9F/8PsNbO5afBN/GT0T6ykjSui2aPuDCwQkeXA58DrqvpWfGOZ1mrqOd3ZUVbFK8t3uJpDVfn9ot8zs2gmNw6+kbsL7raSNq456R8TVXUTMKIZshjD1wZ3YVDXTB54bz2XD+9Gkr/5z6SoL+kZRTOYNmgaPy34qZW0cZWdT2Q8xecT7rq4P1v3VjJ3cfMfq1ZV7iu8r6Gkf3b2z6ykjeusqI3nTBjUiZH5WTz49/XUBJvvSuWqyv8W/i/Prn6WGwbeYCVtPMOK2niOiHD31waws7ya5z5vnkt1qSr3L76fZ1Y/w5SBU7hnzD1W0sYzrKiNJ43r25FzenXgz+9voKo2vnvVqsofl/yRp1c9zXUDruPnY35uJW08xYraeJKIcNfXBlBaUcOzn26J27+jqjyw5AGeWvkU1w24jnvPuddK2niOFbXxrDG9OjC+fw6PfriRiuq6mG9fVXlw6YM8ufJJJvefzC/O+YWVtPEkK2rjaXdd3J/9lXU89fGWmG5XVfm/pf/H9BXTubb/tdx77r02qJLxLHtlGk8bkZ/FxYM785ePNlFWWRuTbdaX9F9W/IVr+l3DL8/9pZW08TR7dRrP+8nF/TlYG+Qv8zed8bbqhyqtL+lfjf2VlbTxPHuFGs8b1DWTy4d346mPt7DnYM0Zbevh5Q/z2BePcXW/q62kTcKwV6lJCD+a2I/quhCPfrDxtLfxyLJHeHT5o1zV9yr+bey/WUmbhGGvVJMQ+uRkcPWoPP762Va+LK8+5ec/svwRHl7+MFf2uZJff+XXVtImodir1SSMOyf0IxRW/vz++lN63mPLH+PhZQ9zRZ8rrKRNQrJXrEkY+R3Sue7sfOYs2h71VWDW7lvLn5f9mct7X85vvvIb/D5vXzjXmKZYUZuE8s8X9UVEePDv0e1VzyyaSVogjXvG3GMlbRKWFbVJKF3bpTHtnB7MXVLMptKDJ1x3X/U+Xt/0Ot/o/Q3apbRrpoTGxF7URS0ifhFZKiKvxTOQMSfzvQv6kBLw88B7J96rnrtuLrXhWm4YdEMzJTMmPk5lj/pOoCheQYyJVk7bFG4Z15NXv9jJmi8PNLlOXbiO59Y+x9iuY+mT1aeZExoTW1EVtYjkAZcB0+Mbx5jo/L/xvclIDvDHd9c1ufy9re9RUlnCtMHTmjmZMbEX7R71A8DPgPDxVhCRO0SkUEQKS0tLYxLOmOPJSk/m9vN78/aq3awoLj9m+cyimXRv253zcs9zIZ0xsXXSohaRy4ESVV18ovVU9XFVLVDVgpycnJgFNOZ4vn1eT7LSk7jvnbVHzF+5ZyXLS5dzw6Ab7Jxp0yJE8yoeB1whIluA54CLRGRGXFMZE4W2qUl896t9+HBdKYu27GuYP7NoJm2S2nBlnytdTGdM7Jy0qFX156qap6o9geuBf6iqHfgznnDT2B5kZ6Rw39trUVVKK0t5a8tbXNX3KjKSM9yOZ0xM2O+FJqGlJwf4/oV9WLh5H59s3MsL614gFA4xZeAUt6MZEzOnVNSq+oGqXh6vMMacjiljutO1XSp/eGclc9bO4fy88+mR2cPtWMbEjO1Rm4SXmuTnBxf1Y2XZfPZV72PqwKluRzImpqyoTYvwrdG5ZHT6lECoC+d0OdftOMbElBW1aRFW71tBKGk7B0vO5Z3Vu92OY0xMWVGbFmFG0QzaJrclP+k87n93HaGwuh3JmJixojYJ78tDX/Le1ve4pt813HXxcNaXHOSV5TvcjmVMzFhRm4Q3Z+0cFOX6gddzydAuDOqayQPvracudNwRD4xJKFbUJqFVB6t5cd2LXJh/IbkZufh8wl0X92fr3krmLi52O54xMWFFbRLaG5vfoKymjKmDDp+SN2FQJ0bkZ/Hg39dTEwy5mM6Y2LCiNglLVZlRNIP+7ftT0LmgYb6IcPfX+rOzvJrnPt/uYkJjYsOK2iSswt2FrN+/nqmDpiIiRyw7r282Y3p14M/vb6Cq1vaqTWKzojYJa8bqGWSlZHFpr0uPWebsVQ+gtKKGv362pfnDGRNDVtQmIRVXFPNB8Qd8q/+3SA2kNrnOmF4dGN8/h0c+2MjBmmAzJzQmdqyoTUJ6bs1zCMJ1A6474Xp3Xdyf/ZV1TP3LZ6ze2fT1FY3xOitqk3Aq6yqZt34eE3tMpEubLidcd0R+Fg/dMIri/VVc8ecF/M9ba6ius2PWJrFYUZuE8+rGV6moq2DaoOiuX3HZ8K6895OvctVZuTzywUYmPfARn2zYE+eUxsSOFbVJKGENM3PNTIZ0HMKInBFRP699m2Tuu3YEM28/BwVumL6Qn76wnLLK2viFNSZGorm4baqIfC4iy0VklYj8ujmCGdOUz3Z+xubyzU2ekheNcX2zeftH4/nuV/swb+kOJt7/Ia8s34mqDeJkvCuaPeoa4CJVHQGMBCaJiA34a1wxo2gGHVM78vWeXz/tbaQm+bnnkoG88s/j6JaVxg9nL+W2ZwrZUVYVw6TGxE40F7dVVT0YeZgUudnuh2l2Ww9sZf6O+UweMJlkf/IZb29It3b87Z/G8cvLBvHpxr1cfP+HPLlgsw2RajwnqmPUIuIXkWVACfCuqi5sYp07RKRQRApLS0tjndMYZhXNIuALMHnA5Jht0+8Tbj+/N+/8eDxn9+zAb15bzdWPfMKaL+1UPuMdURW1qoZUdSSQB4wRkaFNrPO4qhaoakFOTk6sc5pW7mDtQV7a8BKTek4iOy075tvP75DO07eezZ+uH0nxvkouf3ABv7dT+YxHnOpVyMuAD4BJcUljzHG8tOElKoOVR4ySF2siwpUjcxtO5Xu4/lS+jXYqn3FXNGd95IhIVmQ6DZgIrIl3MGPqhTXMrDWzGJEzgqHZx/wyF3PHnMr3l4X87EU7lc+4J5o96q7A+yLyBbAI5xj1a/GNZcxh84vns71ie9RvcImVcX2zeetO51S+uUucU/letVP5jAskHi+6goICLSwsjPl2Tev0nXe+w6byTbx1zVsk+ZJcybBqZzk/n7eCL4rL6dspg+F57RiW69wGd8skPTngSi7TcojIYlUtaGqZvbqMp20s28hnuz7jh2f90LWSBudUvnnf+wqzP9/GP9aU8NG6Pcxb4lxA1yfQJyeDYVbeJk7slWQ8bWbRTJJ9yVzT/xq3oxDw+7hxbE9uHNsTVWX3gRpW7ChnxY5yVu4oZ/76Jso7tx1Dc9sxLK8dg7tm0ibFvuXMqbNXjfGs8ppyXt34Kpf1vowOqR3cjnMEEaFLu1S6tEvl4sGdG+bvPlDNiuLD5b1gwx7mLd0ReQ706JBOhzbJZKUnk5WWRLv0JLLSkmnfJol2aUkN87Mi89umBvD5Tv2t8qZlsaI2njVv/TyqQ9VxPSUv1jpnptJ5cCoTG5V3yYHqhj3v9SUHKa+so6SimnW7KyivrKPiBBc18AlkpiXRPj2ZzLQk0pJ8JAf8JPuF5ICPJL+PZL+PpIBznxy5T4pMJ/mFlMh6Pp/gF8HnA59Iw83vc37w1C9rmK5/jCACgrPs8DQQWeZMRZZHlkmjZUdrPN95RhPzT/DzqfFzTrTt5lL/T/p8Qp+cjJhv34raeFIwHGT2mtkUdC5gQIcBbsc5I50yU5mQmcqEQZ2bXF4XClNeVUdZZR3lVbWUVTrTZVV1lFXWHjFdU+esWxcMUxsKUxcKUxt07muChx/bu+DdkZ2RQuEvJ8Z8u1bUxpM+2P4Buw7t4l/O/he3o8Rdkt9HdkYK2RkpMdtmKKzURsq8vshDYUUVQqqEVQmHlbBCWPWYZapKKOwsC6uC8x/OpEbuaThVUSP/a1imRw4I1PjssiPnN06tx5l/pBP9DHLjzEltlCgl4I/Lv2FFbTxpZtFMurXpxgX5F7gdJSH5fUJasp804lMcpnnZhQOM56zZt4bC3YVMGTgFv8+KxhgrauM5s4pmkRZI45v9vul2FGM8wYraeMq+6n28vul1vtH7G7RLaed2HGM8wYraeMrza5+nNlzLDYNucDuKMZ5hRW08o7KukplFMxmfN54+WX3cjmOMZ1hRG8+Yu34uZTVlfGfYd9yOYoynWFEbT6gL1fHMqmcY3Xk0IzuNdDuOMZ5iRW084bVNr7G7cje3D7vd7SjGeI4VtXFdKBziiZVPMKjDIMZ1G+d2HGM8J5pLceWLyPsiUiQiq0TkzuYIZlqP97a9x9YDW7lt2G2IGyPqGONx0byFPAjcpapLRKQtsFhE3lXV1XHOZloBVeWJFU/QM7MnE7vHfjAbY1qCk+5Rq+ouVV0Sma4AioDceAczrcPHOz+maF8R3x76bXu7uDHHcUrHqEWkJ3AWsLCJZXeISKGIFJaWlsYmnWnxpq+YTuf0zlze+3K3oxjjWVEXtYhkAHOBH6nqgaOXq+rjqlqgqgU5OTmxzGhaqKUlS1m8ezE3D7mZJL9710M0xuuiKmoRScIp6ZmqOi++kUxrMX3FdLJSsrimn/vXQzTGy6I560OAJ4AiVb0//pFMa7B231o+Kv6IqYOmkp6U7nYcYzwtmj3qccCNwEUisixyuzTOuUwL98TKJ0gPpDNl4BS3oxjjeSc9PU9VF8AJriRpzCnadmAbb295m5sH32xDmRoTBXtnoml2T616ioAEuHHwjW5HMSYhWFGbZlVSWcLLG17myr5XkpNuZwcZEw0ratOsnl31LCENcevQW92OYkzCsKI2zaa8ppzn1z3PpJ6TyG+b73YcYxKGFbVpNrOKZlEVrOK2Ybe5HcWYhGJFbZpFZV0lM9fM5IK8C+jfvr/bcYxJKFbUplm8uO5FymvKbW/amNNgRW3irjZUyzOrnqGgc4FdZsuY02BFbeLu1Y2vUlJVYhetNeY0WVGbuAqFQzy58kkGdRjE2G5j3Y5jTEKyojZx9e7Wd9lWsY3vDP+OXWbLmNNkRW3iRlWZvmI6PTN7MqH7BLfjGJOwrKhN3CzYsYC1+9fy7aHfxif2UjPmdNl3j4mb6Sum06VNF7vMljFnyIraxMWS3UtYUrKEW4bcYpfZMuYMWVGbuJi+YjrtU9pzdb+r3Y5iTMKL5lJcT4pIiYisbI5AJvGt2beG+TvmM23wNNICaW7HMSbhRbNH/TQwKc45TAvyxIonaJPUhusHXu92FGNahJMWtap+BOxrhiymBdh6YCvvbH2HyQMmk5mc6XYcY1qEmB2jFpE7RKRQRApLS0tjtVmTYJ5a6Vxm66bBN7kdxZgWI2ZFraqPq2qBqhbk5Nglllqj3Yd28/LGl/lmv2+SnZbtdhxjWgw768PEzLOrn0VVuWXILW5HMaZFsaI2MVFWXcYL617gkl6XkNc2z+04xrQo0ZyeNxv4FBggIsUiYiO/m2PMWhO5zNZQe3kYE2uBk62gqlOaI4hJXIfqDjGzaCYX5F9A3/Z9z2xjB0uhZBV0HwuBlNgENCbBnbSojTmZF9e9yIHaA9w+7PbT30jFbvjkQVj0BASrIDULhl0LZ02FriPBhkg1rZgVtTkj9ZfZGtNlDCNyRpz6Bip2w8d/gsInIVQDwybDgEug6BVY8iws+gt0GuIU9vDroI2dTWJaHytqc0Ze3vgypVWl/Md5/3FqTzywyynoxU9BqM4p4fF3Q8c+zvIhV0HVflg5F5bOhLd/Ae/+CvpPgpFTod/FYIM9mVbCitqctmA4yFMrn2JIxyGM7RrlZbYO7IQFD8DipyEchBFTYPxd0KH3seumtYezb3duu1fDspnwxRxY8xq0yXHK/axp0GlQTD8uY7zGitqclrLqMn772W/ZXrGdP17wx5NfZqt8Byz4Iyx5BjTsFPT5d0GHXtH9g50Hw9d/BxP/Hda/65T2wkfh0z9Dt1Ew8gYY9i2n3I1pYURVY77RgoICLSwsjPl2jTd8sP0D/v2Tf6e8tpzvj/w+tw297fhFXV4M8++HpX91CnrkVKeg2/c48yCH9sAXzzulvXsl+FNg4GXO8ezeF4LPf+b/hjHNREQWq2pBU8tsj9pEraK2gt8v+j0vbXiJ/u3789jFjzGgw4CmVy7bFinoGc7js6bB+T+BrO6xCxB2p3IAAAqtSURBVNQmG8b+E5z7Pdi13CnsFS/AqnmQmQsjrnd+MNQf9zYmQdketYnKZ7s+41cf/4rdlbu5behtfG/E95q+csv+rTD/f2HZLOeUurNuhPN+DFn5zRM0WANr33D+/Q3vOXvx+ec6h0aGXAWp7ZonhzGn6ER71FbU5oQq6yp5YMkDzF4zm56ZPfndeb9jeM7wY1fcv6VRQftg1M1w3o+gnYtvJz+wC754zjlrZO96CKQ6h0ZGTHEOjfjtF0rjHVbU5rQsK1nGvQvuZVvFNqYNmsYPR/3w2Cu27NsM8++D5c+B+GH0LU5BZ3ZzJXOTVGHHElg+G1a+6Jz2l9HZeUPNiCnQZajbCY2xojanpjZUy0PLHuLpVU/TJb0Lvx33W8Z0HXN4hWANrHvL2Xte/65zPvPoW2HcnZDZ1b3g0QjWwPp3nB8s695yThHsMswp7GHXQkYntxOaVsqK2kStaG8Rv1jwCzaUbeCaftdwd8HdZCRnOHulO5c65Vy/V9q2q/MHu3O+C227uB391B3a67yhZvks52MTP/Sd6HxMAy6FpFS3E5pWxM76MCcVDAeZvmI6jy1/jPap7XlowkOMzxsPFV86428smw2lRZHjvJfDyCmJfwpcm45wzh3OrWSNczx7+RxY/zaktIOh33T2tPPPsbFGjKtsj9qwqWwT9y64l5V7V3JJr0u4d9RdtNv6yVFnTpzjlNaQb0JaltuR4yccgs0fOYdGil6Bukpo38v52EdcB+17up3QtFB26MM0qTZUy+w1s3lwyYOkJ6Xzy3438PVd651DG9Xlh89FHnEDZJ/h8KWJqKYCil51fmBtme/M6zEO+k6Ajv0gu5/z1ncbjtXEgBV1K1cTqmFz+WY2lm1suG0q38S2im2ENcwF6fn8W0kJ2aXrIZAGg77hnHfca3xiH9qIpbJtzrsgv5gDe9Ydni8+5008HftGyrvv4enMbnbIxETtjItaRCYBfwL8wHRV/e8TrW9F7Y6qYBVbyrewoWwDm8o3sXH/OjbuW09x5ZeEcb7OfoTukkyfEPSpqWZ4+R7Or6pCuo91ynnwVZCa6fJH4nHVB2DvBti70Tk/e8/6w4/rDh1eLyndeVdkx35OeWf3cw6jpLZzPscpmZCUZmVugDP8Y6KI+IGHgIuBYmCRiLyiqqtjG7P1CmuYmlANNXXVVNcdpLr2INW1B6iurXQe1x2kuq6S6rpKauqqqA5WUh2spjpURVVdFTuqdrOhcjc7Qoeo/7EbUKVHXR0Da+u4rC5I77o6+tbW0UNSSGrb1TmNrlM3GNQXhl5tb7M+FamZkDvKuTWm6owOuHdDpMA3ONM7l8Dql5xj/UfzBSClrVPa9eV9xHTbI+cHUsCf7LxZx58MvqTDjxtPNyyL3HxJzm9H4rMfDAkomrM+xgAbVHUTgIg8B1wJxLyor3tqJNVNvZhPw4l+T4jmYI9G/q9HPG50L0duSAXCODdFI/cQapgXuZem5p3eN05AlVRVugSDDK0NcoUk0yc5i77pXcnP7EFSZq7z63fbrofvbW85fkSgXa5z6/3VI5cFa5w3B5VtdfbIa8qdY+DVB6DmwJHT5cVQEpmuPgAaikNW3wlucuRjJFLu0qjkG82r/9iRyMOj5h39OTr84ATLmlh+wnWjfF68pXeEb78Z881GU9S5wPZGj4uBc45eSUTuAO4A6N799Abe6ZXcntpwMGY/8OUEX7ATf5nlqMdHzq3PJw1LDj/2ieBD8IkPHyAi+PEdXiaCEJkn9esJKb4kUv0ppAZSnHt/KqmBVFICaaQF0kgNpJGS1IbUpDakJqWRkpRBUlK6s7eU3sF5p50NpO9dgRToNNC5nQpV58yT+iIP1UCoFkJBCNcdng7VRh7X32qdN/M0TNc529LwCW4nWq407JnUTzfsqGijeXrUvCM+mCM/ruMta3L5CdaN+nnNIE47QtEUdVOddsxnQ1UfBx4H5xj16YT576nvn87TjGm5RCC5jXNLxDcVmZjwRbFOMdB46LM8YGd84hhjjDlaNEW9COgnIr1EJBm4HnglvrGMMcbUO+mhD1UNisg/A2/jnJ73pKquinsyY4wxQJRjfajqG8Abcc5ijDGmCdEc+jDGGOMiK2pjjPE4K2pjjPE4K2pjjPG4uIyeJyKlwNbTfHo2sCeGcWLN6/nAMsaC1/OB9zN6PR94K2MPVc1pakFcivpMiEjh8UaQ8gKv5wPLGAtezwfez+j1fJAYGcEOfRhjjOdZURtjjMd5sagfdzvASXg9H1jGWPB6PvB+Rq/ng8TI6L1j1MYYY47kxT1qY4wxjVhRG2OMx3mmqEVkkoisFZENInKP23kARCRfRN4XkSIRWSUid0bmdxCRd0VkfeS+vcs5/SKyVERe82i+LBF5UUTWRD6XY72UUUR+HPn6rhSR2SKS6nY+EXlSREpEZGWjecfNJCI/j3zvrBWRr7uY8Q+Rr/MXIvI3EclyK2NT+Rotu1tEVESy3cp3KjxR1I0uoHsJMBiYIiKD3U0FQBC4S1UHAecC34/kugf4u6r2A/4eeeymO4GiRo+9lu9PwFuqOhAYgZPVExlFJBf4IVCgqkNxhvK93gP5ngYmHTWvyUyR1+T1wJDIcx6OfE+5kfFdYKiqDgfWAT93MWNT+RCRfJyLdW9rNM+tz2FUPFHUNLqArqrWAvUX0HWVqu5S1SWR6QqcgsnFyfZMZLVngKvcSQgikgdcBkxvNNtL+TKB8cATAKpaq6pleCgjznC/aSISANJxrmDkaj5V/QjYd9Ts42W6EnhOVWtUdTOwAed7qtkzquo7qhqMPPwM54pQrmQ8zucQ4I/AzzjykoKufA6j5ZWibuoCurkuZWmSiPQEzgIWAp1VdRc4ZQ50ci8ZD+C86Bpfvt1L+XoDpcBTkcMz00WkjVcyquoO4D6cvatdQLmqvuOVfEc5Xiavfv98G6i/JLcnMorIFcAOVV1+1CJP5DserxR1VBfQdYuIZABzgR+p6gG389QTkcuBElVd7HaWEwgAo4BHVPUs4BDuH4ppEDnOeyXQC+gGtBGRae6mOmWe+/4RkXtxDh3OrJ/VxGrNmlFE0oF7gV81tbiJeZ7pIK8UtWcvoCsiSTglPVNV50Vm7xaRrpHlXYESl+KNA64QkS04h4suEpEZHsoHzte2WFUXRh6/iFPcXsk4EdisqqWqWgfMA77ioXyNHS+Tp75/RORm4HJgqh5+o4YXMvbB+YG8PPI9kwcsEZEuHsl3XF4pak9eQFdEBOfYapGq3t9o0SvAzZHpm4GXmzsbgKr+XFXzVLUnzufsH6o6zSv5AFT1S2C7iAyIzJoArMY7GbcB54pIeuTrPQHnbxFeydfY8TK9AlwvIiki0gvoB3zuQj5EZBLwL8AVqlrZaJHrGVV1hap2UtWeke+ZYmBU5DXqer4TUlVP3IBLcf5KvBG41+08kUzn4fz68wWwLHK7FOiI81f39ZH7Dh7IegHwWmTaU/mAkUBh5PP4EtDeSxmBXwNrgJXAX4EUt/MBs3GOmdfhFMptJ8qE8yv9RmAtcImLGTfgHOut/3551K2MTeU7avkWINvNz2G0N3sLuTHGeJxXDn0YY4w5DitqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxuP8PvpAo84fDzDwAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig, ax = plt.subplots()\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][0, :])\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][1, :])\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][2, :])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Split intregrate specification from model-user based specification" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Create middleware function `create_fun` which translates user defined integration function with vars/pars as individual arguments to a scipy compatible input" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 15, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def create_fun(initial_states, parameters):\n", | |
| " \"\"\"Convert integrate statement to scipy-compatible function\"\"\"\n", | |
| " \n", | |
| " # OPTION TO ADD CHECKS on integrate function arguments\n", | |
| " # versus the defined integrate function\n", | |
| " # ... (not done for the example)\n", | |
| "\n", | |
| " def func(t, y, *pars):\n", | |
| " #states = dict(zip(initial_states, y))\n", | |
| " return integrate(t, *y, **parameters) \n", | |
| "\n", | |
| " return func" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 78, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "# User specifies...\n", | |
| "\n", | |
| "# ... state transitions\n", | |
| "def integrate(t, S, I, R, beta, gamma): # All variables and parameters... will be long list or arguments(!)\n", | |
| " \"\"\"Basic SIR model\"\"\"\n", | |
| " N = S + I + R\n", | |
| " dS = -beta*S*I/N\n", | |
| " dI = beta*S*I/N - gamma*I\n", | |
| " dR = gamma*I\n", | |
| "\n", | |
| " return dS, dI, dR\n", | |
| "\n", | |
| "# ... time, parameters and initial conditions\n", | |
| "time = [0, 150]\n", | |
| "parameters = {\"beta\": 0.5, \"gamma\": 0.3}\n", | |
| "initial_states = {\"S\": 7900000, \"I\": 10, \"R\": 0}\n", | |
| "\n", | |
| "# -> prepares the model function and runs with scipy\n", | |
| "fun = create_fun(initial_states, parameters)\n", | |
| "output = solve_ivp(fun, time, list(initial_states.values()), args=list(parameters.values()))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 18, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "[<matplotlib.lines.Line2D at 0x7f4d76649790>]" | |
| ] | |
| }, | |
| "execution_count": 18, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3/8ddnZrISQoCELQn7vgsRpShVwRaXqtWKIrhV66+9va1ttb229vbetrd3ab3WeutaXMsiKtR9b13ABQmbLGFfA0jCkhDIOjOf3x9nEgIEGGAm50zyeT4c58w5Zw7vJJN3Tk7OfI+oKsYYY7zL53YAY4wxJ2ZFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHhe3ohaRJ0WkRERWRrn+ZBFZLSKrRGRWvHIZY0yikXidRy0i44GDwLOqOvQk6/YDngcuUtX9ItJJVUviEswYYxJM3PaoVfUjYF/jeSLSR0TeEpHFIjJfRAZGFn0HeEhV90eeayVtjDERzX2M+nHgB6o6GrgbeDgyvz/QX0Q+FpHPRGRSM+cyxhjPCjTXPyQiGcBXgBdEpH52SqMc/YALgDxgvogMVdWy5spnjDFe1WxFjbP3XqaqI5tYVgx8pqp1wGYRWYtT3IuaMZ8xxnhSsx36UNUDOCV8LYA4RkQWvwRcGJmfjXMoZFNzZTPGGC+L5+l5s4FPgQEiUiwitwFTgdtEZDmwCrgysvrbwF4RWQ28D/xUVffGK5sxxiSSuJ2eZ4wxJjbsnYnGGONxcfljYnZ2tvbs2TMemzbGmBZp8eLFe1Q1p6llcSnqnj17UlhYGI9NG2NMiyQiW4+3zA59GGOMx1lRG2OMx1lRG2OMx0VV1CLy48jwoytFZLaIpMY7mDHGGMdJi1pEcoEfAgWR4Ur9wPXxDmaMMcYR7aGPAJAmIgEgHdgZv0jGGGMaO2lRq+oO4D5gG7ALKFfVd45eT0TuEJFCESksLS2NfVJjjGmlTnoetYi0xxmToxdQhjNM6TRVndF4PVV9HGe8aQoKCk7rfekP/n09wVD4ZIFOZ9OHnx7F5gRBxFm3fnn90KzOfDk8v9E8n08I+AR/5L7xY3/DtA+/D/w+n7OOCEl+ISM1QGZqEm1TA2SkBJAz/DiNMS1HNG94mQhsVtVSABGZhzOu9IwTPus0PPrhRqrqQsdd3lqGJfEJtE1NIjPNKe/MxtNpxz4e2KUt+R3S3Y5tjImTaIp6G3CuiKQDVcAEIC5vO1z9m+a9sEtTA1KpgkaW6RHztOEHRcN9ZF79+uEwBMNhQqqEwkowpIRVCYadx/W34BHTYYIh5WBNkIrqOg5UBTlQXceBqjoOVAcj93Vs2VPZMP9Q7bE/zHp2TOe8ftmc3y+HsX06kpmaFJfPmTGm+Z20qFV1oYi8CCwBgsBSIoc4El1ThxcOz/LuoYdgKExFtVPo+yvrWLptP/PX72Hekh3M+Gwbfp8wIq8d5/fL4fx+2YzIzyLJb6fMG5Oo4jLMaUFBgdpYH82vNhhmybb9LFi/h/kb9rCiuIywQkZKgLF9OnJ+v2zO65tNr+w2dgzcGI8RkcWqWtDkMivqlqusspZPNu5l/vo9zF9fSvH+KgBys9Kc0o4Ud1Z6sstJjTFW1AZVZeveSuZv2MP8daV8unEvFTVB2qYE+Ovt5zAyP8vtiMa0albU5hjBUJhl28v4yfPLKa+q47k7zmVQ10y3YxnTap2oqO0vTK1UwO+joGcHZt5+DmlJfm584nM2lR50O5YxpglW1K1cfod0Ztx+DqrKtOkLKd5f6XYkY8xRrKgNfTtl8OxtY6ioCTJt+kJKKqrdjmSMacSK2gAwpFs7nr51DCUVNdw4/XP2H6p1O5IxJsKK2jQY3aM9028qYPPeQ9z81OdUVNe5HckYgxW1OcpX+mbz8A2jWL3zALc9U0hVE29XN8Y0Lytqc4yJgztz/3UjWbRlH9+dsZiaoJW1MW6yojZNumJEN/776mF8uK6UO2cvO/nws8aYuLGiNsd13dnd+dfLB/PWqi/52dwvCIdbyTizxnhMNMOcmlbstvN6cagmyP3vrqNNcoDfXDnEBnQypplZUZuT+sFFfTlUE+SxjzbRJiXAv0waYGVtTDOyojYnJSLcc8lADtYEefTDjbRNDfD9C/u6HcuYVsOK2kRFRPjtlUOprA3xh7fXkp7s59ZxvdyOZUyrEM3FbQcAcxrN6g38SlUfiFsq40k+n/CHbw2nsjbIr19dTZuUAJML8t2OZUyLd9KzPlR1raqOVNWRwGigEvhb3JMZTwr4fTw45SzO75fNPXO/4PUvdrkdyZgW71RPz5sAbFTVrfEIYxJDSsDP4zcWMLpHe+58bikbbXhUY+LqVIv6emB2UwtE5A4RKRSRwtLS0jNPZjwtLdnPQ1NHAfD8ou0upzGmZYu6qEUkGbgCeKGp5ar6uKoWqGpBTk5OrPIZD+vUNpWLBnZi7pId1Nk7F42Jm1PZo74EWKKqu+MVxiSeyQX57DlYw/trStyOYkyLdSpFPYXjHPYwrdcFA3LIaZvC84XFbkcxpsWKqqhFJB24GJgX3zgm0QT8Pq4elcv7a0vsyjDGxElURa2qlaraUVXL4x3IJJ7JBfmEwsrfluxwO4oxLZKNnmfOWJ+cDAp6tGdO4XZUbYQ9Y2LNitrExOSCfDaVHmLJtv1uRzGmxbGiNjFx6fCupCf7eX6R/VHRmFizojYxkZES4PLhXXnti50cqgm6HceYFsWK2sTM5IJ8DtWGeH2Fjf9hTCxZUZuYGd2jPb2z2/BCob2l3JhYsqI2MSMiXFuQz6It+9lkAzUZEzNW1CamrhmVi98nvLDY/qhoTKxYUZuY6pSZyoUDcpi7uJigDdRkTExYUZuYu7Ygn5KKGj5cZ8PdGhMLVtQm5i4a2InsjGSetz8qGhMTVtQm5pL8Pq4elcffi0rYc7DG7TjGJDwrahMX147OI2gDNRkTE1bUJi76dW7LWd2zeN4GajLmjFlRm7iZXJDP+pKDLNte5nYUYxKaFbWJm8uHdyUtyW9XfzHmDEV7hZcsEXlRRNaISJGIjI13MJP42qYmcemwrry6fCeVtTZQkzGnK9o96j8Bb6nqQGAEUBS/SKYlmVyQx8GaIG+u+NLtKMYkrJMWtYhkAuOBJwBUtVZV7aCjicqYXh3o2THdzqk25gxEs0fdGygFnhKRpSIyXUTaHL2SiNwhIoUiUlhaau9IM476gZoWbt7Hlj2H3I5jTEKKpqgDwCjgEVU9CzgE3HP0Sqr6uKoWqGpBTk5OjGOaRHbNqDx8Ai8str1qY05HNEVdDBSr6sLI4xdxituYqHRpl8pX++fw4uJiQmE7p9qYU3XSolbVL4HtIjIgMmsCsDquqUyLM7kgn90HavhovR0WM+ZURXvWxw+AmSLyBTAS+M/4RTIt0YRBnenQJtmu/mLMaQhEs5KqLgMK4pzFtGDJAR/fPCuXZz/dwt6DNXTMSHE7kjEJw96ZaJrN5IJ86kLKS8t2uh3FmIRiRW2azYAubRmR144XbKAmY06JFbVpVpPPzmfNlxWs2FHudhRjEoYVtWlW3xjRjZSAjzmL7I+KxkTLito0q8zIQE2vLNtJVW3I7TjGJAQratPsri3Io6ImyNurbKAmY6JhRW2a3bm9OtK9gw3UZEy0rKhNs/P5hGtH5/HJxr1s21vpdhxjPM+K2rjimtF5iMCLNlCTMScV1TsTjYm1bllpnN/PGajpzon98fvE7UjmKKpKSEOENdxwA5xpwqgqqnp4GiWsh6cb7htNOxumYbp+WcN0o3WOyHL0jEi+Y+Y1sV5UH+tpPu9oAQnQO6t3TLZ1xHZjvkVjonRdQT7fn7WEjzfsYXx/Gxr3ZA7WHuRA7QGqglVU1lVSGaw8fN9ouqqu6shldZVUBauoC9cR0hDBcJBgONgwHQqHCGrT88yp6ZjakQ+u+yDm27WiNq6ZOLgTWelJzCncbkUNVAWr2HlwJzsO7nBuFTsaposPFlNRWxHVdlL8KaQH0klPSictkNZwn+HLIOALEJAAfp+fgC+AX/wk+ZLwi79hXkACzjKf35kvfnziQ0Tw4dwL0jCvYRpx1mk0LTi/KdVPi0Qec/g3qCPWO846kRnHOGad48yLRv2/eyaSfclnvI2mWFEb16QE/Fw1MpdZC7ex/1At7dvE50XuFapK8cFiiiuKGwq5+GBxQynvrd57xPop/hS6ZXQjNyOX4TnDyc3IJSsli7SkNKeII2Xc+D41kErAZ9/WLY19RY2rJhfk8/QnW3h52Q5uGdfL7ThxsalsE69vfp03N7/J9orDfzz1i58ubbqQl5HHV/O/Sm5G7hG37LTsmOzlmcRnRW1cNbhbJsNy2/F8YXGLKupdB3fx5pY3eWPTG6zdvxaf+BjTZQw3D76Z3lm9yc3IpVN6J9v7NVGxV4lx3eSCPP715VWs3FHO0Nx2bsc5bfur9/POlnd4Y/MbLClZAsDw7OHcM+Yevt7z62SnZbuc0CSqqIpaRLYAFUAICKqqXUTAxMwVI3L57etFPF+4PeGK+lDdIf6x7R+8uflNPt35KUEN0rtdb35w1g+4pOcl5Gfmux3RtACnskd9oaruiVsS02q1S09i0pAuvLR0B7+4dBCpSX63I51QbaiWBTsW8ObmN/lg+wdUh6rp2qYrNw25iUt7XUr/9v3t2LKJKTv0YTzhurPzeWX5Tt5ZvZsrRnRzO06TVJWHlj3ErDWzqKitoH1Ke67seyWX9rqUkZ1G4hN7o6+Jj2iLWoF3RESBx1T18aNXEJE7gDsAunfvHruEplUY27sjuVlpPL9ouyeLWlX53cLfMWftHCZ2n8jV/a7m3G7nkuRLcjuaaQWi3QUYp6qjgEuA74vI+KNXUNXHVbVAVQtycuzNC+bU+HzCtQV5fLxxD9v3eWugJlXlPxf+J3PWzuHWIbdy/wX3c37e+VbSptlEVdSqujNyXwL8DRgTz1CmdfrW6DwA5i4pdjnJYarKf33+Xzy39jluGXILPx79Yzv+bJrdSYtaRNqISNv6aeBrwMp4BzOtT177dM7rm80LhcWEwu5f/FZV+Z9F/8PsNbO5afBN/GT0T6ykjSui2aPuDCwQkeXA58DrqvpWfGOZ1mrqOd3ZUVbFK8t3uJpDVfn9ot8zs2gmNw6+kbsL7raSNq456R8TVXUTMKIZshjD1wZ3YVDXTB54bz2XD+9Gkr/5z6SoL+kZRTOYNmgaPy34qZW0cZWdT2Q8xecT7rq4P1v3VjJ3cfMfq1ZV7iu8r6Gkf3b2z6ykjeusqI3nTBjUiZH5WTz49/XUBJvvSuWqyv8W/i/Prn6WGwbeYCVtPMOK2niOiHD31waws7ya5z5vnkt1qSr3L76fZ1Y/w5SBU7hnzD1W0sYzrKiNJ43r25FzenXgz+9voKo2vnvVqsofl/yRp1c9zXUDruPnY35uJW08xYraeJKIcNfXBlBaUcOzn26J27+jqjyw5AGeWvkU1w24jnvPuddK2niOFbXxrDG9OjC+fw6PfriRiuq6mG9fVXlw6YM8ufJJJvefzC/O+YWVtPEkK2rjaXdd3J/9lXU89fGWmG5XVfm/pf/H9BXTubb/tdx77r02qJLxLHtlGk8bkZ/FxYM785ePNlFWWRuTbdaX9F9W/IVr+l3DL8/9pZW08TR7dRrP+8nF/TlYG+Qv8zed8bbqhyqtL+lfjf2VlbTxPHuFGs8b1DWTy4d346mPt7DnYM0Zbevh5Q/z2BePcXW/q62kTcKwV6lJCD+a2I/quhCPfrDxtLfxyLJHeHT5o1zV9yr+bey/WUmbhGGvVJMQ+uRkcPWoPP762Va+LK8+5ec/svwRHl7+MFf2uZJff+XXVtImodir1SSMOyf0IxRW/vz++lN63mPLH+PhZQ9zRZ8rrKRNQrJXrEkY+R3Sue7sfOYs2h71VWDW7lvLn5f9mct7X85vvvIb/D5vXzjXmKZYUZuE8s8X9UVEePDv0e1VzyyaSVogjXvG3GMlbRKWFbVJKF3bpTHtnB7MXVLMptKDJ1x3X/U+Xt/0Ot/o/Q3apbRrpoTGxF7URS0ifhFZKiKvxTOQMSfzvQv6kBLw88B7J96rnrtuLrXhWm4YdEMzJTMmPk5lj/pOoCheQYyJVk7bFG4Z15NXv9jJmi8PNLlOXbiO59Y+x9iuY+mT1aeZExoTW1EVtYjkAZcB0+Mbx5jo/L/xvclIDvDHd9c1ufy9re9RUlnCtMHTmjmZMbEX7R71A8DPgPDxVhCRO0SkUEQKS0tLYxLOmOPJSk/m9vN78/aq3awoLj9m+cyimXRv253zcs9zIZ0xsXXSohaRy4ESVV18ovVU9XFVLVDVgpycnJgFNOZ4vn1eT7LSk7jvnbVHzF+5ZyXLS5dzw6Ab7Jxp0yJE8yoeB1whIluA54CLRGRGXFMZE4W2qUl896t9+HBdKYu27GuYP7NoJm2S2nBlnytdTGdM7Jy0qFX156qap6o9geuBf6iqHfgznnDT2B5kZ6Rw39trUVVKK0t5a8tbXNX3KjKSM9yOZ0xM2O+FJqGlJwf4/oV9WLh5H59s3MsL614gFA4xZeAUt6MZEzOnVNSq+oGqXh6vMMacjiljutO1XSp/eGclc9bO4fy88+mR2cPtWMbEjO1Rm4SXmuTnBxf1Y2XZfPZV72PqwKluRzImpqyoTYvwrdG5ZHT6lECoC+d0OdftOMbElBW1aRFW71tBKGk7B0vO5Z3Vu92OY0xMWVGbFmFG0QzaJrclP+k87n93HaGwuh3JmJixojYJ78tDX/Le1ve4pt813HXxcNaXHOSV5TvcjmVMzFhRm4Q3Z+0cFOX6gddzydAuDOqayQPvracudNwRD4xJKFbUJqFVB6t5cd2LXJh/IbkZufh8wl0X92fr3krmLi52O54xMWFFbRLaG5vfoKymjKmDDp+SN2FQJ0bkZ/Hg39dTEwy5mM6Y2LCiNglLVZlRNIP+7ftT0LmgYb6IcPfX+rOzvJrnPt/uYkJjYsOK2iSswt2FrN+/nqmDpiIiRyw7r282Y3p14M/vb6Cq1vaqTWKzojYJa8bqGWSlZHFpr0uPWebsVQ+gtKKGv362pfnDGRNDVtQmIRVXFPNB8Qd8q/+3SA2kNrnOmF4dGN8/h0c+2MjBmmAzJzQmdqyoTUJ6bs1zCMJ1A6474Xp3Xdyf/ZV1TP3LZ6ze2fT1FY3xOitqk3Aq6yqZt34eE3tMpEubLidcd0R+Fg/dMIri/VVc8ecF/M9ba6ius2PWJrFYUZuE8+rGV6moq2DaoOiuX3HZ8K6895OvctVZuTzywUYmPfARn2zYE+eUxsSOFbVJKGENM3PNTIZ0HMKInBFRP699m2Tuu3YEM28/BwVumL6Qn76wnLLK2viFNSZGorm4baqIfC4iy0VklYj8ujmCGdOUz3Z+xubyzU2ekheNcX2zeftH4/nuV/swb+kOJt7/Ia8s34mqDeJkvCuaPeoa4CJVHQGMBCaJiA34a1wxo2gGHVM78vWeXz/tbaQm+bnnkoG88s/j6JaVxg9nL+W2ZwrZUVYVw6TGxE40F7dVVT0YeZgUudnuh2l2Ww9sZf6O+UweMJlkf/IZb29It3b87Z/G8cvLBvHpxr1cfP+HPLlgsw2RajwnqmPUIuIXkWVACfCuqi5sYp07RKRQRApLS0tjndMYZhXNIuALMHnA5Jht0+8Tbj+/N+/8eDxn9+zAb15bzdWPfMKaL+1UPuMdURW1qoZUdSSQB4wRkaFNrPO4qhaoakFOTk6sc5pW7mDtQV7a8BKTek4iOy075tvP75DO07eezZ+uH0nxvkouf3ABv7dT+YxHnOpVyMuAD4BJcUljzHG8tOElKoOVR4ySF2siwpUjcxtO5Xu4/lS+jXYqn3FXNGd95IhIVmQ6DZgIrIl3MGPqhTXMrDWzGJEzgqHZx/wyF3PHnMr3l4X87EU7lc+4J5o96q7A+yLyBbAI5xj1a/GNZcxh84vns71ie9RvcImVcX2zeetO51S+uUucU/letVP5jAskHi+6goICLSwsjPl2Tev0nXe+w6byTbx1zVsk+ZJcybBqZzk/n7eCL4rL6dspg+F57RiW69wGd8skPTngSi7TcojIYlUtaGqZvbqMp20s28hnuz7jh2f90LWSBudUvnnf+wqzP9/GP9aU8NG6Pcxb4lxA1yfQJyeDYVbeJk7slWQ8bWbRTJJ9yVzT/xq3oxDw+7hxbE9uHNsTVWX3gRpW7ChnxY5yVu4oZ/76Jso7tx1Dc9sxLK8dg7tm0ibFvuXMqbNXjfGs8ppyXt34Kpf1vowOqR3cjnMEEaFLu1S6tEvl4sGdG+bvPlDNiuLD5b1gwx7mLd0ReQ706JBOhzbJZKUnk5WWRLv0JLLSkmnfJol2aUkN87Mi89umBvD5Tv2t8qZlsaI2njVv/TyqQ9VxPSUv1jpnptJ5cCoTG5V3yYHqhj3v9SUHKa+so6SimnW7KyivrKPiBBc18AlkpiXRPj2ZzLQk0pJ8JAf8JPuF5ICPJL+PZL+PpIBznxy5T4pMJ/mFlMh6Pp/gF8HnA59Iw83vc37w1C9rmK5/jCACgrPs8DQQWeZMRZZHlkmjZUdrPN95RhPzT/DzqfFzTrTt5lL/T/p8Qp+cjJhv34raeFIwHGT2mtkUdC5gQIcBbsc5I50yU5mQmcqEQZ2bXF4XClNeVUdZZR3lVbWUVTrTZVV1lFXWHjFdU+esWxcMUxsKUxcKUxt07muChx/bu+DdkZ2RQuEvJ8Z8u1bUxpM+2P4Buw7t4l/O/he3o8Rdkt9HdkYK2RkpMdtmKKzURsq8vshDYUUVQqqEVQmHlbBCWPWYZapKKOwsC6uC8x/OpEbuaThVUSP/a1imRw4I1PjssiPnN06tx5l/pBP9DHLjzEltlCgl4I/Lv2FFbTxpZtFMurXpxgX5F7gdJSH5fUJasp804lMcpnnZhQOM56zZt4bC3YVMGTgFv8+KxhgrauM5s4pmkRZI45v9vul2FGM8wYraeMq+6n28vul1vtH7G7RLaed2HGM8wYraeMrza5+nNlzLDYNucDuKMZ5hRW08o7KukplFMxmfN54+WX3cjmOMZ1hRG8+Yu34uZTVlfGfYd9yOYoynWFEbT6gL1fHMqmcY3Xk0IzuNdDuOMZ5iRW084bVNr7G7cje3D7vd7SjGeI4VtXFdKBziiZVPMKjDIMZ1G+d2HGM8J5pLceWLyPsiUiQiq0TkzuYIZlqP97a9x9YDW7lt2G2IGyPqGONx0byFPAjcpapLRKQtsFhE3lXV1XHOZloBVeWJFU/QM7MnE7vHfjAbY1qCk+5Rq+ouVV0Sma4AioDceAczrcPHOz+maF8R3x76bXu7uDHHcUrHqEWkJ3AWsLCJZXeISKGIFJaWlsYmnWnxpq+YTuf0zlze+3K3oxjjWVEXtYhkAHOBH6nqgaOXq+rjqlqgqgU5OTmxzGhaqKUlS1m8ezE3D7mZJL9710M0xuuiKmoRScIp6ZmqOi++kUxrMX3FdLJSsrimn/vXQzTGy6I560OAJ4AiVb0//pFMa7B231o+Kv6IqYOmkp6U7nYcYzwtmj3qccCNwEUisixyuzTOuUwL98TKJ0gPpDNl4BS3oxjjeSc9PU9VF8AJriRpzCnadmAbb295m5sH32xDmRoTBXtnoml2T616ioAEuHHwjW5HMSYhWFGbZlVSWcLLG17myr5XkpNuZwcZEw0ratOsnl31LCENcevQW92OYkzCsKI2zaa8ppzn1z3PpJ6TyG+b73YcYxKGFbVpNrOKZlEVrOK2Ybe5HcWYhGJFbZpFZV0lM9fM5IK8C+jfvr/bcYxJKFbUplm8uO5FymvKbW/amNNgRW3irjZUyzOrnqGgc4FdZsuY02BFbeLu1Y2vUlJVYhetNeY0WVGbuAqFQzy58kkGdRjE2G5j3Y5jTEKyojZx9e7Wd9lWsY3vDP+OXWbLmNNkRW3iRlWZvmI6PTN7MqH7BLfjGJOwrKhN3CzYsYC1+9fy7aHfxif2UjPmdNl3j4mb6Sum06VNF7vMljFnyIraxMWS3UtYUrKEW4bcYpfZMuYMWVGbuJi+YjrtU9pzdb+r3Y5iTMKL5lJcT4pIiYisbI5AJvGt2beG+TvmM23wNNICaW7HMSbhRbNH/TQwKc45TAvyxIonaJPUhusHXu92FGNahJMWtap+BOxrhiymBdh6YCvvbH2HyQMmk5mc6XYcY1qEmB2jFpE7RKRQRApLS0tjtVmTYJ5a6Vxm66bBN7kdxZgWI2ZFraqPq2qBqhbk5Nglllqj3Yd28/LGl/lmv2+SnZbtdhxjWgw768PEzLOrn0VVuWXILW5HMaZFsaI2MVFWXcYL617gkl6XkNc2z+04xrQo0ZyeNxv4FBggIsUiYiO/m2PMWhO5zNZQe3kYE2uBk62gqlOaI4hJXIfqDjGzaCYX5F9A3/Z9z2xjB0uhZBV0HwuBlNgENCbBnbSojTmZF9e9yIHaA9w+7PbT30jFbvjkQVj0BASrIDULhl0LZ02FriPBhkg1rZgVtTkj9ZfZGtNlDCNyRpz6Bip2w8d/gsInIVQDwybDgEug6BVY8iws+gt0GuIU9vDroI2dTWJaHytqc0Ze3vgypVWl/Md5/3FqTzywyynoxU9BqM4p4fF3Q8c+zvIhV0HVflg5F5bOhLd/Ae/+CvpPgpFTod/FYIM9mVbCitqctmA4yFMrn2JIxyGM7RrlZbYO7IQFD8DipyEchBFTYPxd0KH3seumtYezb3duu1fDspnwxRxY8xq0yXHK/axp0GlQTD8uY7zGitqclrLqMn772W/ZXrGdP17wx5NfZqt8Byz4Iyx5BjTsFPT5d0GHXtH9g50Hw9d/BxP/Hda/65T2wkfh0z9Dt1Ew8gYY9i2n3I1pYURVY77RgoICLSwsjPl2jTd8sP0D/v2Tf6e8tpzvj/w+tw297fhFXV4M8++HpX91CnrkVKeg2/c48yCH9sAXzzulvXsl+FNg4GXO8ezeF4LPf+b/hjHNREQWq2pBU8tsj9pEraK2gt8v+j0vbXiJ/u3789jFjzGgw4CmVy7bFinoGc7js6bB+T+BrO6xCxB2p3IAAAqtSURBVNQmG8b+E5z7Pdi13CnsFS/AqnmQmQsjrnd+MNQf9zYmQdketYnKZ7s+41cf/4rdlbu5behtfG/E95q+csv+rTD/f2HZLOeUurNuhPN+DFn5zRM0WANr33D+/Q3vOXvx+ec6h0aGXAWp7ZonhzGn6ER71FbU5oQq6yp5YMkDzF4zm56ZPfndeb9jeM7wY1fcv6VRQftg1M1w3o+gnYtvJz+wC754zjlrZO96CKQ6h0ZGTHEOjfjtF0rjHVbU5rQsK1nGvQvuZVvFNqYNmsYPR/3w2Cu27NsM8++D5c+B+GH0LU5BZ3ZzJXOTVGHHElg+G1a+6Jz2l9HZeUPNiCnQZajbCY2xojanpjZUy0PLHuLpVU/TJb0Lvx33W8Z0HXN4hWANrHvL2Xte/65zPvPoW2HcnZDZ1b3g0QjWwPp3nB8s695yThHsMswp7GHXQkYntxOaVsqK2kStaG8Rv1jwCzaUbeCaftdwd8HdZCRnOHulO5c65Vy/V9q2q/MHu3O+C227uB391B3a67yhZvks52MTP/Sd6HxMAy6FpFS3E5pWxM76MCcVDAeZvmI6jy1/jPap7XlowkOMzxsPFV86428smw2lRZHjvJfDyCmJfwpcm45wzh3OrWSNczx7+RxY/zaktIOh33T2tPPPsbFGjKtsj9qwqWwT9y64l5V7V3JJr0u4d9RdtNv6yVFnTpzjlNaQb0JaltuR4yccgs0fOYdGil6Bukpo38v52EdcB+17up3QtFB26MM0qTZUy+w1s3lwyYOkJ6Xzy3438PVd651DG9Xlh89FHnEDZJ/h8KWJqKYCil51fmBtme/M6zEO+k6Ajv0gu5/z1ncbjtXEgBV1K1cTqmFz+WY2lm1suG0q38S2im2ENcwF6fn8W0kJ2aXrIZAGg77hnHfca3xiH9qIpbJtzrsgv5gDe9Ydni8+5008HftGyrvv4enMbnbIxETtjItaRCYBfwL8wHRV/e8TrW9F7Y6qYBVbyrewoWwDm8o3sXH/OjbuW09x5ZeEcb7OfoTukkyfEPSpqWZ4+R7Or6pCuo91ynnwVZCa6fJH4nHVB2DvBti70Tk/e8/6w4/rDh1eLyndeVdkx35OeWf3cw6jpLZzPscpmZCUZmVugDP8Y6KI+IGHgIuBYmCRiLyiqqtjG7P1CmuYmlANNXXVVNcdpLr2INW1B6iurXQe1x2kuq6S6rpKauqqqA5WUh2spjpURVVdFTuqdrOhcjc7Qoeo/7EbUKVHXR0Da+u4rC5I77o6+tbW0UNSSGrb1TmNrlM3GNQXhl5tb7M+FamZkDvKuTWm6owOuHdDpMA3ONM7l8Dql5xj/UfzBSClrVPa9eV9xHTbI+cHUsCf7LxZx58MvqTDjxtPNyyL3HxJzm9H4rMfDAkomrM+xgAbVHUTgIg8B1wJxLyor3tqJNVNvZhPw4l+T4jmYI9G/q9HPG50L0duSAXCODdFI/cQapgXuZem5p3eN05AlVRVugSDDK0NcoUk0yc5i77pXcnP7EFSZq7z63fbrofvbW85fkSgXa5z6/3VI5cFa5w3B5VtdfbIa8qdY+DVB6DmwJHT5cVQEpmuPgAaikNW3wlucuRjJFLu0qjkG82r/9iRyMOj5h39OTr84ATLmlh+wnWjfF68pXeEb78Z881GU9S5wPZGj4uBc45eSUTuAO4A6N799Abe6ZXcntpwMGY/8OUEX7ATf5nlqMdHzq3PJw1LDj/2ieBD8IkPHyAi+PEdXiaCEJkn9esJKb4kUv0ppAZSnHt/KqmBVFICaaQF0kgNpJGS1IbUpDakJqWRkpRBUlK6s7eU3sF5p50NpO9dgRToNNC5nQpV58yT+iIP1UCoFkJBCNcdng7VRh7X32qdN/M0TNc529LwCW4nWq407JnUTzfsqGijeXrUvCM+mCM/ruMta3L5CdaN+nnNIE47QtEUdVOddsxnQ1UfBx4H5xj16YT576nvn87TjGm5RCC5jXNLxDcVmZjwRbFOMdB46LM8YGd84hhjjDlaNEW9COgnIr1EJBm4HnglvrGMMcbUO+mhD1UNisg/A2/jnJ73pKquinsyY4wxQJRjfajqG8Abcc5ijDGmCdEc+jDGGOMiK2pjjPE4K2pjjPE4K2pjjPG4uIyeJyKlwNbTfHo2sCeGcWLN6/nAMsaC1/OB9zN6PR94K2MPVc1pakFcivpMiEjh8UaQ8gKv5wPLGAtezwfez+j1fJAYGcEOfRhjjOdZURtjjMd5sagfdzvASXg9H1jGWPB6PvB+Rq/ng8TI6L1j1MYYY47kxT1qY4wxjVhRG2OMx3mmqEVkkoisFZENInKP23kARCRfRN4XkSIRWSUid0bmdxCRd0VkfeS+vcs5/SKyVERe82i+LBF5UUTWRD6XY72UUUR+HPn6rhSR2SKS6nY+EXlSREpEZGWjecfNJCI/j3zvrBWRr7uY8Q+Rr/MXIvI3EclyK2NT+Rotu1tEVESy3cp3KjxR1I0uoHsJMBiYIiKD3U0FQBC4S1UHAecC34/kugf4u6r2A/4eeeymO4GiRo+9lu9PwFuqOhAYgZPVExlFJBf4IVCgqkNxhvK93gP5ngYmHTWvyUyR1+T1wJDIcx6OfE+5kfFdYKiqDgfWAT93MWNT+RCRfJyLdW9rNM+tz2FUPFHUNLqArqrWAvUX0HWVqu5S1SWR6QqcgsnFyfZMZLVngKvcSQgikgdcBkxvNNtL+TKB8cATAKpaq6pleCgjznC/aSISANJxrmDkaj5V/QjYd9Ts42W6EnhOVWtUdTOwAed7qtkzquo7qhqMPPwM54pQrmQ8zucQ4I/AzzjykoKufA6j5ZWibuoCurkuZWmSiPQEzgIWAp1VdRc4ZQ50ci8ZD+C86Bpfvt1L+XoDpcBTkcMz00WkjVcyquoO4D6cvatdQLmqvuOVfEc5Xiavfv98G6i/JLcnMorIFcAOVV1+1CJP5DserxR1VBfQdYuIZABzgR+p6gG389QTkcuBElVd7HaWEwgAo4BHVPUs4BDuH4ppEDnOeyXQC+gGtBGRae6mOmWe+/4RkXtxDh3OrJ/VxGrNmlFE0oF7gV81tbiJeZ7pIK8UtWcvoCsiSTglPVNV50Vm7xaRrpHlXYESl+KNA64QkS04h4suEpEZHsoHzte2WFUXRh6/iFPcXsk4EdisqqWqWgfMA77ioXyNHS+Tp75/RORm4HJgqh5+o4YXMvbB+YG8PPI9kwcsEZEuHsl3XF4pak9eQFdEBOfYapGq3t9o0SvAzZHpm4GXmzsbgKr+XFXzVLUnzufsH6o6zSv5AFT1S2C7iAyIzJoArMY7GbcB54pIeuTrPQHnbxFeydfY8TK9AlwvIiki0gvoB3zuQj5EZBLwL8AVqlrZaJHrGVV1hap2UtWeke+ZYmBU5DXqer4TUlVP3IBLcf5KvBG41+08kUzn4fz68wWwLHK7FOiI81f39ZH7Dh7IegHwWmTaU/mAkUBh5PP4EtDeSxmBXwNrgJXAX4EUt/MBs3GOmdfhFMptJ8qE8yv9RmAtcImLGTfgHOut/3551K2MTeU7avkWINvNz2G0N3sLuTHGeJxXDn0YY4w5DitqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxuP8PvpAo84fDzDwAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig, ax = plt.subplots()\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][0, :])\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][1, :])\n", | |
| "ax.plot(output[\"t\"], output[\"y\"][2, :])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Abstract away the boilerplate in Base class and define a new model as subclass" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "class BaseModel:\n", | |
| " \n", | |
| " def __init__(self, states, parameters):\n", | |
| " \"\"\"\"\"\"\n", | |
| " self.parameters = parameters\n", | |
| " self.initial_states = states\n", | |
| " \n", | |
| " # OPTION TO ADD CHECKS on incoming arguments on initialization\n", | |
| " # versus the defined integrate variables and class attributes of the subclass\n", | |
| " # ... @Joris -> we should add POC here as well, as it needs to compare the \n", | |
| " # incoming dicts with the class attributes of the subclass?\n", | |
| " # (and I rather not add a init to the subclass)\n", | |
| "\n", | |
| " def integrate(self):\n", | |
| " \"\"\"to overwrite in subclasses\"\"\"\n", | |
| " raise NotImplementedError \n", | |
| " \n", | |
| " def create_fun(self):\n", | |
| " \"\"\"Convert integrate statement to scipy-compatible function\"\"\"\n", | |
| "\n", | |
| " def func(t, y, *pars):\n", | |
| " #states = dict(zip(initial_states, y)) # @Joris, added value? \n", | |
| " return self.integrate(t, *y, **self.parameters) \n", | |
| "\n", | |
| " return func\n", | |
| " \n", | |
| " def sim(self, time):\n", | |
| " \"\"\"\"\"\" \n", | |
| " fun = self.create_fun()\n", | |
| " output = solve_ivp(fun, time, \n", | |
| " list(self.initial_states.values()), \n", | |
| " args=list(self.parameters.values()))\n", | |
| " return output[\"t\"], self.array_to_variables(output[\"y\"]) # map back to variable names\n", | |
| " \n", | |
| " def array_to_variables(self, y):\n", | |
| " \"\"\"Convert array (used by scipy) to dictionary (used by model API)\"\"\"\n", | |
| " return dict(zip(self.state_names, y))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "# User (model developer) specifies...\n", | |
| "\n", | |
| "# ... state transitions in a subclass\n", | |
| "class SIR(BaseModel):\n", | |
| "\n", | |
| " # state variables and parameters\n", | |
| " state_names = ['S', 'I', 'R']\n", | |
| " parameter_names = ['beta', 'gamma']\n", | |
| " \n", | |
| " @staticmethod\n", | |
| " def integrate(t, S, I, R, beta, gamma): # All variables and parameters... will be long list or arguments(!)\n", | |
| " \"\"\"Basic SIR model\"\"\"\n", | |
| " N = S + I + R\n", | |
| " dS = -beta*S*I/N\n", | |
| " dI = beta*S*I/N - gamma*I\n", | |
| " dR = gamma*I\n", | |
| " \n", | |
| " return dS, dI, dR\n", | |
| "\n", | |
| "# ... parameters and initial conditions\n", | |
| "parameters = {\"beta\": 0.5, \"gamma\": 0.3}\n", | |
| "initial_states = {\"S\": 7900000, \"I\": 10, \"R\": 0}\n", | |
| "\n", | |
| "# -> user initiates the model\n", | |
| "sir_model = SIR(initial_states, parameters)\n", | |
| "\n", | |
| "# -> user runs a simulation for a defined time period\n", | |
| "time = [0, 150]\n", | |
| "t, output = sir_model.sim(time)\n", | |
| "\n", | |
| "# -> user can do fit, mc,... using the model class instance `sir_model`" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "[<matplotlib.lines.Line2D at 0x7f880e210110>]" | |
| ] | |
| }, | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3/8ddnZrISQoCELQn7vgsRpShVwRaXqtWKIrhV66+9va1ttb229vbetrd3ab3WeutaXMsiKtR9b13ABQmbLGFfA0jCkhDIOjOf3x9nEgIEGGAm50zyeT4c58w5Zw7vJJN3Tk7OfI+oKsYYY7zL53YAY4wxJ2ZFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHhe3ohaRJ0WkRERWRrn+ZBFZLSKrRGRWvHIZY0yikXidRy0i44GDwLOqOvQk6/YDngcuUtX9ItJJVUviEswYYxJM3PaoVfUjYF/jeSLSR0TeEpHFIjJfRAZGFn0HeEhV90eeayVtjDERzX2M+nHgB6o6GrgbeDgyvz/QX0Q+FpHPRGRSM+cyxhjPCjTXPyQiGcBXgBdEpH52SqMc/YALgDxgvogMVdWy5spnjDFe1WxFjbP3XqaqI5tYVgx8pqp1wGYRWYtT3IuaMZ8xxnhSsx36UNUDOCV8LYA4RkQWvwRcGJmfjXMoZFNzZTPGGC+L5+l5s4FPgQEiUiwitwFTgdtEZDmwCrgysvrbwF4RWQ28D/xUVffGK5sxxiSSuJ2eZ4wxJjbsnYnGGONxcfljYnZ2tvbs2TMemzbGmBZp8eLFe1Q1p6llcSnqnj17UlhYGI9NG2NMiyQiW4+3zA59GGOMx1lRG2OMx1lRG2OMx0VV1CLy48jwoytFZLaIpMY7mDHGGMdJi1pEcoEfAgWR4Ur9wPXxDmaMMcYR7aGPAJAmIgEgHdgZv0jGGGMaO2lRq+oO4D5gG7ALKFfVd45eT0TuEJFCESksLS2NfVJjjGmlTnoetYi0xxmToxdQhjNM6TRVndF4PVV9HGe8aQoKCk7rfekP/n09wVD4ZIFOZ9OHnx7F5gRBxFm3fnn90KzOfDk8v9E8n08I+AR/5L7xY3/DtA+/D/w+n7OOCEl+ISM1QGZqEm1TA2SkBJAz/DiNMS1HNG94mQhsVtVSABGZhzOu9IwTPus0PPrhRqrqQsdd3lqGJfEJtE1NIjPNKe/MxtNpxz4e2KUt+R3S3Y5tjImTaIp6G3CuiKQDVcAEIC5vO1z9m+a9sEtTA1KpgkaW6RHztOEHRcN9ZF79+uEwBMNhQqqEwkowpIRVCYadx/W34BHTYYIh5WBNkIrqOg5UBTlQXceBqjoOVAcj93Vs2VPZMP9Q7bE/zHp2TOe8ftmc3y+HsX06kpmaFJfPmTGm+Z20qFV1oYi8CCwBgsBSIoc4El1ThxcOz/LuoYdgKExFtVPo+yvrWLptP/PX72Hekh3M+Gwbfp8wIq8d5/fL4fx+2YzIzyLJb6fMG5Oo4jLMaUFBgdpYH82vNhhmybb9LFi/h/kb9rCiuIywQkZKgLF9OnJ+v2zO65tNr+w2dgzcGI8RkcWqWtDkMivqlqusspZPNu5l/vo9zF9fSvH+KgBys9Kc0o4Ud1Z6sstJjTFW1AZVZeveSuZv2MP8daV8unEvFTVB2qYE+Ovt5zAyP8vtiMa0albU5hjBUJhl28v4yfPLKa+q47k7zmVQ10y3YxnTap2oqO0vTK1UwO+joGcHZt5+DmlJfm584nM2lR50O5YxpglW1K1cfod0Ztx+DqrKtOkLKd5f6XYkY8xRrKgNfTtl8OxtY6ioCTJt+kJKKqrdjmSMacSK2gAwpFs7nr51DCUVNdw4/XP2H6p1O5IxJsKK2jQY3aM9028qYPPeQ9z81OdUVNe5HckYgxW1OcpX+mbz8A2jWL3zALc9U0hVE29XN8Y0Lytqc4yJgztz/3UjWbRlH9+dsZiaoJW1MW6yojZNumJEN/776mF8uK6UO2cvO/nws8aYuLGiNsd13dnd+dfLB/PWqi/52dwvCIdbyTizxnhMNMOcmlbstvN6cagmyP3vrqNNcoDfXDnEBnQypplZUZuT+sFFfTlUE+SxjzbRJiXAv0waYGVtTDOyojYnJSLcc8lADtYEefTDjbRNDfD9C/u6HcuYVsOK2kRFRPjtlUOprA3xh7fXkp7s59ZxvdyOZUyrEM3FbQcAcxrN6g38SlUfiFsq40k+n/CHbw2nsjbIr19dTZuUAJML8t2OZUyLd9KzPlR1raqOVNWRwGigEvhb3JMZTwr4fTw45SzO75fNPXO/4PUvdrkdyZgW71RPz5sAbFTVrfEIYxJDSsDP4zcWMLpHe+58bikbbXhUY+LqVIv6emB2UwtE5A4RKRSRwtLS0jNPZjwtLdnPQ1NHAfD8ou0upzGmZYu6qEUkGbgCeKGp5ar6uKoWqGpBTk5OrPIZD+vUNpWLBnZi7pId1Nk7F42Jm1PZo74EWKKqu+MVxiSeyQX57DlYw/trStyOYkyLdSpFPYXjHPYwrdcFA3LIaZvC84XFbkcxpsWKqqhFJB24GJgX3zgm0QT8Pq4elcv7a0vsyjDGxElURa2qlaraUVXL4x3IJJ7JBfmEwsrfluxwO4oxLZKNnmfOWJ+cDAp6tGdO4XZUbYQ9Y2LNitrExOSCfDaVHmLJtv1uRzGmxbGiNjFx6fCupCf7eX6R/VHRmFizojYxkZES4PLhXXnti50cqgm6HceYFsWK2sTM5IJ8DtWGeH2Fjf9hTCxZUZuYGd2jPb2z2/BCob2l3JhYsqI2MSMiXFuQz6It+9lkAzUZEzNW1CamrhmVi98nvLDY/qhoTKxYUZuY6pSZyoUDcpi7uJigDdRkTExYUZuYu7Ygn5KKGj5cZ8PdGhMLVtQm5i4a2InsjGSetz8qGhMTVtQm5pL8Pq4elcffi0rYc7DG7TjGJDwrahMX147OI2gDNRkTE1bUJi76dW7LWd2zeN4GajLmjFlRm7iZXJDP+pKDLNte5nYUYxKaFbWJm8uHdyUtyW9XfzHmDEV7hZcsEXlRRNaISJGIjI13MJP42qYmcemwrry6fCeVtTZQkzGnK9o96j8Bb6nqQGAEUBS/SKYlmVyQx8GaIG+u+NLtKMYkrJMWtYhkAuOBJwBUtVZV7aCjicqYXh3o2THdzqk25gxEs0fdGygFnhKRpSIyXUTaHL2SiNwhIoUiUlhaau9IM476gZoWbt7Hlj2H3I5jTEKKpqgDwCjgEVU9CzgE3HP0Sqr6uKoWqGpBTk5OjGOaRHbNqDx8Ai8str1qY05HNEVdDBSr6sLI4xdxituYqHRpl8pX++fw4uJiQmE7p9qYU3XSolbVL4HtIjIgMmsCsDquqUyLM7kgn90HavhovR0WM+ZURXvWxw+AmSLyBTAS+M/4RTIt0YRBnenQJtmu/mLMaQhEs5KqLgMK4pzFtGDJAR/fPCuXZz/dwt6DNXTMSHE7kjEJw96ZaJrN5IJ86kLKS8t2uh3FmIRiRW2azYAubRmR144XbKAmY06JFbVpVpPPzmfNlxWs2FHudhRjEoYVtWlW3xjRjZSAjzmL7I+KxkTLito0q8zIQE2vLNtJVW3I7TjGJAQratPsri3Io6ImyNurbKAmY6JhRW2a3bm9OtK9gw3UZEy0rKhNs/P5hGtH5/HJxr1s21vpdhxjPM+K2rjimtF5iMCLNlCTMScV1TsTjYm1bllpnN/PGajpzon98fvE7UjmKKpKSEOENdxwA5xpwqgqqnp4GiWsh6cb7htNOxumYbp+WcN0o3WOyHL0jEi+Y+Y1sV5UH+tpPu9oAQnQO6t3TLZ1xHZjvkVjonRdQT7fn7WEjzfsYXx/Gxr3ZA7WHuRA7QGqglVU1lVSGaw8fN9ouqqu6shldZVUBauoC9cR0hDBcJBgONgwHQqHCGrT88yp6ZjakQ+u+yDm27WiNq6ZOLgTWelJzCncbkUNVAWr2HlwJzsO7nBuFTsaposPFlNRWxHVdlL8KaQH0klPSictkNZwn+HLIOALEJAAfp+fgC+AX/wk+ZLwi79hXkACzjKf35kvfnziQ0Tw4dwL0jCvYRpx1mk0LTi/KdVPi0Qec/g3qCPWO846kRnHOGad48yLRv2/eyaSfclnvI2mWFEb16QE/Fw1MpdZC7ex/1At7dvE50XuFapK8cFiiiuKGwq5+GBxQynvrd57xPop/hS6ZXQjNyOX4TnDyc3IJSsli7SkNKeII2Xc+D41kErAZ9/WLY19RY2rJhfk8/QnW3h52Q5uGdfL7ThxsalsE69vfp03N7/J9orDfzz1i58ubbqQl5HHV/O/Sm5G7hG37LTsmOzlmcRnRW1cNbhbJsNy2/F8YXGLKupdB3fx5pY3eWPTG6zdvxaf+BjTZQw3D76Z3lm9yc3IpVN6J9v7NVGxV4lx3eSCPP715VWs3FHO0Nx2bsc5bfur9/POlnd4Y/MbLClZAsDw7OHcM+Yevt7z62SnZbuc0CSqqIpaRLYAFUAICKqqXUTAxMwVI3L57etFPF+4PeGK+lDdIf6x7R+8uflNPt35KUEN0rtdb35w1g+4pOcl5Gfmux3RtACnskd9oaruiVsS02q1S09i0pAuvLR0B7+4dBCpSX63I51QbaiWBTsW8ObmN/lg+wdUh6rp2qYrNw25iUt7XUr/9v3t2LKJKTv0YTzhurPzeWX5Tt5ZvZsrRnRzO06TVJWHlj3ErDWzqKitoH1Ke67seyWX9rqUkZ1G4hN7o6+Jj2iLWoF3RESBx1T18aNXEJE7gDsAunfvHruEplUY27sjuVlpPL9ouyeLWlX53cLfMWftHCZ2n8jV/a7m3G7nkuRLcjuaaQWi3QUYp6qjgEuA74vI+KNXUNXHVbVAVQtycuzNC+bU+HzCtQV5fLxxD9v3eWugJlXlPxf+J3PWzuHWIbdy/wX3c37e+VbSptlEVdSqujNyXwL8DRgTz1CmdfrW6DwA5i4pdjnJYarKf33+Xzy39jluGXILPx79Yzv+bJrdSYtaRNqISNv6aeBrwMp4BzOtT177dM7rm80LhcWEwu5f/FZV+Z9F/8PsNbO5afBN/GT0T6ykjSui2aPuDCwQkeXA58DrqvpWfGOZ1mrqOd3ZUVbFK8t3uJpDVfn9ot8zs2gmNw6+kbsL7raSNq456R8TVXUTMKIZshjD1wZ3YVDXTB54bz2XD+9Gkr/5z6SoL+kZRTOYNmgaPy34qZW0cZWdT2Q8xecT7rq4P1v3VjJ3cfMfq1ZV7iu8r6Gkf3b2z6ykjeusqI3nTBjUiZH5WTz49/XUBJvvSuWqyv8W/i/Prn6WGwbeYCVtPMOK2niOiHD31waws7ya5z5vnkt1qSr3L76fZ1Y/w5SBU7hnzD1W0sYzrKiNJ43r25FzenXgz+9voKo2vnvVqsofl/yRp1c9zXUDruPnY35uJW08xYraeJKIcNfXBlBaUcOzn26J27+jqjyw5AGeWvkU1w24jnvPuddK2niOFbXxrDG9OjC+fw6PfriRiuq6mG9fVXlw6YM8ufJJJvefzC/O+YWVtPEkK2rjaXdd3J/9lXU89fGWmG5XVfm/pf/H9BXTubb/tdx77r02qJLxLHtlGk8bkZ/FxYM785ePNlFWWRuTbdaX9F9W/IVr+l3DL8/9pZW08TR7dRrP+8nF/TlYG+Qv8zed8bbqhyqtL+lfjf2VlbTxPHuFGs8b1DWTy4d346mPt7DnYM0Zbevh5Q/z2BePcXW/q62kTcKwV6lJCD+a2I/quhCPfrDxtLfxyLJHeHT5o1zV9yr+bey/WUmbhGGvVJMQ+uRkcPWoPP762Va+LK8+5ec/svwRHl7+MFf2uZJff+XXVtImodir1SSMOyf0IxRW/vz++lN63mPLH+PhZQ9zRZ8rrKRNQrJXrEkY+R3Sue7sfOYs2h71VWDW7lvLn5f9mct7X85vvvIb/D5vXzjXmKZYUZuE8s8X9UVEePDv0e1VzyyaSVogjXvG3GMlbRKWFbVJKF3bpTHtnB7MXVLMptKDJ1x3X/U+Xt/0Ot/o/Q3apbRrpoTGxF7URS0ifhFZKiKvxTOQMSfzvQv6kBLw88B7J96rnrtuLrXhWm4YdEMzJTMmPk5lj/pOoCheQYyJVk7bFG4Z15NXv9jJmi8PNLlOXbiO59Y+x9iuY+mT1aeZExoTW1EVtYjkAZcB0+Mbx5jo/L/xvclIDvDHd9c1ufy9re9RUlnCtMHTmjmZMbEX7R71A8DPgPDxVhCRO0SkUEQKS0tLYxLOmOPJSk/m9vN78/aq3awoLj9m+cyimXRv253zcs9zIZ0xsXXSohaRy4ESVV18ovVU9XFVLVDVgpycnJgFNOZ4vn1eT7LSk7jvnbVHzF+5ZyXLS5dzw6Ab7Jxp0yJE8yoeB1whIluA54CLRGRGXFMZE4W2qUl896t9+HBdKYu27GuYP7NoJm2S2nBlnytdTGdM7Jy0qFX156qap6o9geuBf6iqHfgznnDT2B5kZ6Rw39trUVVKK0t5a8tbXNX3KjKSM9yOZ0xM2O+FJqGlJwf4/oV9WLh5H59s3MsL614gFA4xZeAUt6MZEzOnVNSq+oGqXh6vMMacjiljutO1XSp/eGclc9bO4fy88+mR2cPtWMbEjO1Rm4SXmuTnBxf1Y2XZfPZV72PqwKluRzImpqyoTYvwrdG5ZHT6lECoC+d0OdftOMbElBW1aRFW71tBKGk7B0vO5Z3Vu92OY0xMWVGbFmFG0QzaJrclP+k87n93HaGwuh3JmJixojYJ78tDX/Le1ve4pt813HXxcNaXHOSV5TvcjmVMzFhRm4Q3Z+0cFOX6gddzydAuDOqayQPvracudNwRD4xJKFbUJqFVB6t5cd2LXJh/IbkZufh8wl0X92fr3krmLi52O54xMWFFbRLaG5vfoKymjKmDDp+SN2FQJ0bkZ/Hg39dTEwy5mM6Y2LCiNglLVZlRNIP+7ftT0LmgYb6IcPfX+rOzvJrnPt/uYkJjYsOK2iSswt2FrN+/nqmDpiIiRyw7r282Y3p14M/vb6Cq1vaqTWKzojYJa8bqGWSlZHFpr0uPWebsVQ+gtKKGv362pfnDGRNDVtQmIRVXFPNB8Qd8q/+3SA2kNrnOmF4dGN8/h0c+2MjBmmAzJzQmdqyoTUJ6bs1zCMJ1A6474Xp3Xdyf/ZV1TP3LZ6ze2fT1FY3xOitqk3Aq6yqZt34eE3tMpEubLidcd0R+Fg/dMIri/VVc8ecF/M9ba6ius2PWJrFYUZuE8+rGV6moq2DaoOiuX3HZ8K6895OvctVZuTzywUYmPfARn2zYE+eUxsSOFbVJKGENM3PNTIZ0HMKInBFRP699m2Tuu3YEM28/BwVumL6Qn76wnLLK2viFNSZGorm4baqIfC4iy0VklYj8ujmCGdOUz3Z+xubyzU2ekheNcX2zeftH4/nuV/swb+kOJt7/Ia8s34mqDeJkvCuaPeoa4CJVHQGMBCaJiA34a1wxo2gGHVM78vWeXz/tbaQm+bnnkoG88s/j6JaVxg9nL+W2ZwrZUVYVw6TGxE40F7dVVT0YeZgUudnuh2l2Ww9sZf6O+UweMJlkf/IZb29It3b87Z/G8cvLBvHpxr1cfP+HPLlgsw2RajwnqmPUIuIXkWVACfCuqi5sYp07RKRQRApLS0tjndMYZhXNIuALMHnA5Jht0+8Tbj+/N+/8eDxn9+zAb15bzdWPfMKaL+1UPuMdURW1qoZUdSSQB4wRkaFNrPO4qhaoakFOTk6sc5pW7mDtQV7a8BKTek4iOy075tvP75DO07eezZ+uH0nxvkouf3ABv7dT+YxHnOpVyMuAD4BJcUljzHG8tOElKoOVR4ySF2siwpUjcxtO5Xu4/lS+jXYqn3FXNGd95IhIVmQ6DZgIrIl3MGPqhTXMrDWzGJEzgqHZx/wyF3PHnMr3l4X87EU7lc+4J5o96q7A+yLyBbAI5xj1a/GNZcxh84vns71ie9RvcImVcX2zeetO51S+uUucU/letVP5jAskHi+6goICLSwsjPl2Tev0nXe+w6byTbx1zVsk+ZJcybBqZzk/n7eCL4rL6dspg+F57RiW69wGd8skPTngSi7TcojIYlUtaGqZvbqMp20s28hnuz7jh2f90LWSBudUvnnf+wqzP9/GP9aU8NG6Pcxb4lxA1yfQJyeDYVbeJk7slWQ8bWbRTJJ9yVzT/xq3oxDw+7hxbE9uHNsTVWX3gRpW7ChnxY5yVu4oZ/76Jso7tx1Dc9sxLK8dg7tm0ibFvuXMqbNXjfGs8ppyXt34Kpf1vowOqR3cjnMEEaFLu1S6tEvl4sGdG+bvPlDNiuLD5b1gwx7mLd0ReQ706JBOhzbJZKUnk5WWRLv0JLLSkmnfJol2aUkN87Mi89umBvD5Tv2t8qZlsaI2njVv/TyqQ9VxPSUv1jpnptJ5cCoTG5V3yYHqhj3v9SUHKa+so6SimnW7KyivrKPiBBc18AlkpiXRPj2ZzLQk0pJ8JAf8JPuF5ICPJL+PZL+PpIBznxy5T4pMJ/mFlMh6Pp/gF8HnA59Iw83vc37w1C9rmK5/jCACgrPs8DQQWeZMRZZHlkmjZUdrPN95RhPzT/DzqfFzTrTt5lL/T/p8Qp+cjJhv34raeFIwHGT2mtkUdC5gQIcBbsc5I50yU5mQmcqEQZ2bXF4XClNeVUdZZR3lVbWUVTrTZVV1lFXWHjFdU+esWxcMUxsKUxcKUxt07muChx/bu+DdkZ2RQuEvJ8Z8u1bUxpM+2P4Buw7t4l/O/he3o8Rdkt9HdkYK2RkpMdtmKKzURsq8vshDYUUVQqqEVQmHlbBCWPWYZapKKOwsC6uC8x/OpEbuaThVUSP/a1imRw4I1PjssiPnN06tx5l/pBP9DHLjzEltlCgl4I/Lv2FFbTxpZtFMurXpxgX5F7gdJSH5fUJasp804lMcpnnZhQOM56zZt4bC3YVMGTgFv8+KxhgrauM5s4pmkRZI45v9vul2FGM8wYraeMq+6n28vul1vtH7G7RLaed2HGM8wYraeMrza5+nNlzLDYNucDuKMZ5hRW08o7KukplFMxmfN54+WX3cjmOMZ1hRG8+Yu34uZTVlfGfYd9yOYoynWFEbT6gL1fHMqmcY3Xk0IzuNdDuOMZ5iRW084bVNr7G7cje3D7vd7SjGeI4VtXFdKBziiZVPMKjDIMZ1G+d2HGM8J5pLceWLyPsiUiQiq0TkzuYIZlqP97a9x9YDW7lt2G2IGyPqGONx0byFPAjcpapLRKQtsFhE3lXV1XHOZloBVeWJFU/QM7MnE7vHfjAbY1qCk+5Rq+ouVV0Sma4AioDceAczrcPHOz+maF8R3x76bXu7uDHHcUrHqEWkJ3AWsLCJZXeISKGIFJaWlsYmnWnxpq+YTuf0zlze+3K3oxjjWVEXtYhkAHOBH6nqgaOXq+rjqlqgqgU5OTmxzGhaqKUlS1m8ezE3D7mZJL9710M0xuuiKmoRScIp6ZmqOi++kUxrMX3FdLJSsrimn/vXQzTGy6I560OAJ4AiVb0//pFMa7B231o+Kv6IqYOmkp6U7nYcYzwtmj3qccCNwEUisixyuzTOuUwL98TKJ0gPpDNl4BS3oxjjeSc9PU9VF8AJriRpzCnadmAbb295m5sH32xDmRoTBXtnoml2T616ioAEuHHwjW5HMSYhWFGbZlVSWcLLG17myr5XkpNuZwcZEw0ratOsnl31LCENcevQW92OYkzCsKI2zaa8ppzn1z3PpJ6TyG+b73YcYxKGFbVpNrOKZlEVrOK2Ybe5HcWYhGJFbZpFZV0lM9fM5IK8C+jfvr/bcYxJKFbUplm8uO5FymvKbW/amNNgRW3irjZUyzOrnqGgc4FdZsuY02BFbeLu1Y2vUlJVYhetNeY0WVGbuAqFQzy58kkGdRjE2G5j3Y5jTEKyojZx9e7Wd9lWsY3vDP+OXWbLmNNkRW3iRlWZvmI6PTN7MqH7BLfjGJOwrKhN3CzYsYC1+9fy7aHfxif2UjPmdNl3j4mb6Sum06VNF7vMljFnyIraxMWS3UtYUrKEW4bcYpfZMuYMWVGbuJi+YjrtU9pzdb+r3Y5iTMKL5lJcT4pIiYisbI5AJvGt2beG+TvmM23wNNICaW7HMSbhRbNH/TQwKc45TAvyxIonaJPUhusHXu92FGNahJMWtap+BOxrhiymBdh6YCvvbH2HyQMmk5mc6XYcY1qEmB2jFpE7RKRQRApLS0tjtVmTYJ5a6Vxm66bBN7kdxZgWI2ZFraqPq2qBqhbk5Nglllqj3Yd28/LGl/lmv2+SnZbtdhxjWgw768PEzLOrn0VVuWXILW5HMaZFsaI2MVFWXcYL617gkl6XkNc2z+04xrQo0ZyeNxv4FBggIsUiYiO/m2PMWhO5zNZQe3kYE2uBk62gqlOaI4hJXIfqDjGzaCYX5F9A3/Z9z2xjB0uhZBV0HwuBlNgENCbBnbSojTmZF9e9yIHaA9w+7PbT30jFbvjkQVj0BASrIDULhl0LZ02FriPBhkg1rZgVtTkj9ZfZGtNlDCNyRpz6Bip2w8d/gsInIVQDwybDgEug6BVY8iws+gt0GuIU9vDroI2dTWJaHytqc0Ze3vgypVWl/Md5/3FqTzywyynoxU9BqM4p4fF3Q8c+zvIhV0HVflg5F5bOhLd/Ae/+CvpPgpFTod/FYIM9mVbCitqctmA4yFMrn2JIxyGM7RrlZbYO7IQFD8DipyEchBFTYPxd0KH3seumtYezb3duu1fDspnwxRxY8xq0yXHK/axp0GlQTD8uY7zGitqclrLqMn772W/ZXrGdP17wx5NfZqt8Byz4Iyx5BjTsFPT5d0GHXtH9g50Hw9d/BxP/Hda/65T2wkfh0z9Dt1Ew8gYY9i2n3I1pYURVY77RgoICLSwsjPl2jTd8sP0D/v2Tf6e8tpzvj/w+tw297fhFXV4M8++HpX91CnrkVKeg2/c48yCH9sAXzzulvXsl+FNg4GXO8ezeF4LPf+b/hjHNREQWq2pBU8tsj9pEraK2gt8v+j0vbXiJ/u3789jFjzGgw4CmVy7bFinoGc7js6bB+T+BrO6xCxB2p3IAAAqtSURBVNQmG8b+E5z7Pdi13CnsFS/AqnmQmQsjrnd+MNQf9zYmQdketYnKZ7s+41cf/4rdlbu5behtfG/E95q+csv+rTD/f2HZLOeUurNuhPN+DFn5zRM0WANr33D+/Q3vOXvx+ec6h0aGXAWp7ZonhzGn6ER71FbU5oQq6yp5YMkDzF4zm56ZPfndeb9jeM7wY1fcv6VRQftg1M1w3o+gnYtvJz+wC754zjlrZO96CKQ6h0ZGTHEOjfjtF0rjHVbU5rQsK1nGvQvuZVvFNqYNmsYPR/3w2Cu27NsM8++D5c+B+GH0LU5BZ3ZzJXOTVGHHElg+G1a+6Jz2l9HZeUPNiCnQZajbCY2xojanpjZUy0PLHuLpVU/TJb0Lvx33W8Z0HXN4hWANrHvL2Xte/65zPvPoW2HcnZDZ1b3g0QjWwPp3nB8s695yThHsMswp7GHXQkYntxOaVsqK2kStaG8Rv1jwCzaUbeCaftdwd8HdZCRnOHulO5c65Vy/V9q2q/MHu3O+C227uB391B3a67yhZvks52MTP/Sd6HxMAy6FpFS3E5pWxM76MCcVDAeZvmI6jy1/jPap7XlowkOMzxsPFV86428smw2lRZHjvJfDyCmJfwpcm45wzh3OrWSNczx7+RxY/zaktIOh33T2tPPPsbFGjKtsj9qwqWwT9y64l5V7V3JJr0u4d9RdtNv6yVFnTpzjlNaQb0JaltuR4yccgs0fOYdGil6Bukpo38v52EdcB+17up3QtFB26MM0qTZUy+w1s3lwyYOkJ6Xzy3438PVd651DG9Xlh89FHnEDZJ/h8KWJqKYCil51fmBtme/M6zEO+k6Ajv0gu5/z1ncbjtXEgBV1K1cTqmFz+WY2lm1suG0q38S2im2ENcwF6fn8W0kJ2aXrIZAGg77hnHfca3xiH9qIpbJtzrsgv5gDe9Ydni8+5008HftGyrvv4enMbnbIxETtjItaRCYBfwL8wHRV/e8TrW9F7Y6qYBVbyrewoWwDm8o3sXH/OjbuW09x5ZeEcb7OfoTukkyfEPSpqWZ4+R7Or6pCuo91ynnwVZCa6fJH4nHVB2DvBti70Tk/e8/6w4/rDh1eLyndeVdkx35OeWf3cw6jpLZzPscpmZCUZmVugDP8Y6KI+IGHgIuBYmCRiLyiqqtjG7P1CmuYmlANNXXVVNcdpLr2INW1B6iurXQe1x2kuq6S6rpKauqqqA5WUh2spjpURVVdFTuqdrOhcjc7Qoeo/7EbUKVHXR0Da+u4rC5I77o6+tbW0UNSSGrb1TmNrlM3GNQXhl5tb7M+FamZkDvKuTWm6owOuHdDpMA3ONM7l8Dql5xj/UfzBSClrVPa9eV9xHTbI+cHUsCf7LxZx58MvqTDjxtPNyyL3HxJzm9H4rMfDAkomrM+xgAbVHUTgIg8B1wJxLyor3tqJNVNvZhPw4l+T4jmYI9G/q9HPG50L0duSAXCODdFI/cQapgXuZem5p3eN05AlVRVugSDDK0NcoUk0yc5i77pXcnP7EFSZq7z63fbrofvbW85fkSgXa5z6/3VI5cFa5w3B5VtdfbIa8qdY+DVB6DmwJHT5cVQEpmuPgAaikNW3wlucuRjJFLu0qjkG82r/9iRyMOj5h39OTr84ATLmlh+wnWjfF68pXeEb78Z881GU9S5wPZGj4uBc45eSUTuAO4A6N799Abe6ZXcntpwMGY/8OUEX7ATf5nlqMdHzq3PJw1LDj/2ieBD8IkPHyAi+PEdXiaCEJkn9esJKb4kUv0ppAZSnHt/KqmBVFICaaQF0kgNpJGS1IbUpDakJqWRkpRBUlK6s7eU3sF5p50NpO9dgRToNNC5nQpV58yT+iIP1UCoFkJBCNcdng7VRh7X32qdN/M0TNc529LwCW4nWq407JnUTzfsqGijeXrUvCM+mCM/ruMta3L5CdaN+nnNIE47QtEUdVOddsxnQ1UfBx4H5xj16YT576nvn87TjGm5RCC5jXNLxDcVmZjwRbFOMdB46LM8YGd84hhjjDlaNEW9COgnIr1EJBm4HnglvrGMMcbUO+mhD1UNisg/A2/jnJ73pKquinsyY4wxQJRjfajqG8Abcc5ijDGmCdEc+jDGGOMiK2pjjPE4K2pjjPE4K2pjjPG4uIyeJyKlwNbTfHo2sCeGcWLN6/nAMsaC1/OB9zN6PR94K2MPVc1pakFcivpMiEjh8UaQ8gKv5wPLGAtezwfez+j1fJAYGcEOfRhjjOdZURtjjMd5sagfdzvASXg9H1jGWPB6PvB+Rq/ng8TI6L1j1MYYY47kxT1qY4wxjVhRG2OMx3mmqEVkkoisFZENInKP23kARCRfRN4XkSIRWSUid0bmdxCRd0VkfeS+vcs5/SKyVERe82i+LBF5UUTWRD6XY72UUUR+HPn6rhSR2SKS6nY+EXlSREpEZGWjecfNJCI/j3zvrBWRr7uY8Q+Rr/MXIvI3EclyK2NT+Rotu1tEVESy3cp3KjxR1I0uoHsJMBiYIiKD3U0FQBC4S1UHAecC34/kugf4u6r2A/4eeeymO4GiRo+9lu9PwFuqOhAYgZPVExlFJBf4IVCgqkNxhvK93gP5ngYmHTWvyUyR1+T1wJDIcx6OfE+5kfFdYKiqDgfWAT93MWNT+RCRfJyLdW9rNM+tz2FUPFHUNLqArqrWAvUX0HWVqu5S1SWR6QqcgsnFyfZMZLVngKvcSQgikgdcBkxvNNtL+TKB8cATAKpaq6pleCgjznC/aSISANJxrmDkaj5V/QjYd9Ts42W6EnhOVWtUdTOwAed7qtkzquo7qhqMPPwM54pQrmQ8zucQ4I/AzzjykoKufA6j5ZWibuoCurkuZWmSiPQEzgIWAp1VdRc4ZQ50ci8ZD+C86Bpfvt1L+XoDpcBTkcMz00WkjVcyquoO4D6cvatdQLmqvuOVfEc5Xiavfv98G6i/JLcnMorIFcAOVV1+1CJP5DserxR1VBfQdYuIZABzgR+p6gG389QTkcuBElVd7HaWEwgAo4BHVPUs4BDuH4ppEDnOeyXQC+gGtBGRae6mOmWe+/4RkXtxDh3OrJ/VxGrNmlFE0oF7gV81tbiJeZ7pIK8UtWcvoCsiSTglPVNV50Vm7xaRrpHlXYESl+KNA64QkS04h4suEpEZHsoHzte2WFUXRh6/iFPcXsk4EdisqqWqWgfMA77ioXyNHS+Tp75/RORm4HJgqh5+o4YXMvbB+YG8PPI9kwcsEZEuHsl3XF4pak9eQFdEBOfYapGq3t9o0SvAzZHpm4GXmzsbgKr+XFXzVLUnzufsH6o6zSv5AFT1S2C7iAyIzJoArMY7GbcB54pIeuTrPQHnbxFeydfY8TK9AlwvIiki0gvoB3zuQj5EZBLwL8AVqlrZaJHrGVV1hap2UtWeke+ZYmBU5DXqer4TUlVP3IBLcf5KvBG41+08kUzn4fz68wWwLHK7FOiI81f39ZH7Dh7IegHwWmTaU/mAkUBh5PP4EtDeSxmBXwNrgJXAX4EUt/MBs3GOmdfhFMptJ8qE8yv9RmAtcImLGTfgHOut/3551K2MTeU7avkWINvNz2G0N3sLuTHGeJxXDn0YY4w5DitqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxuP8PvpAo84fDzDwAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig, ax = plt.subplots()\n", | |
| "ax.plot(t, output[\"S\"])\n", | |
| "ax.plot(t, output[\"I\"])\n", | |
| "ax.plot(t, output[\"R\"])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python [conda env:COVID_MODEL]", | |
| "language": "python", | |
| "name": "conda-env-COVID_MODEL-py" | |
| }, | |
| "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.6" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Beautiful!