diff --git a/modelhub/notebooks/19_05_2020_Africa_SEIRD.ipynb b/modelhub/notebooks/19_05_2020_Africa_SEIRD.ipynb new file mode 100644 index 0000000..d773e18 --- /dev/null +++ b/modelhub/notebooks/19_05_2020_Africa_SEIRD.ipynb @@ -0,0 +1,465 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.integrate import odeint\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline \n", + "\n", + "#!pip install mpld3\n", + "import mpld3\n", + "mpld3.enable_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def plotseird(t, S, E, I, R, D=None, L=None, R0=None, Alpha=None, CFR=None):\n", + " f, ax = plt.subplots(1,1,figsize=(10,4))\n", + " ax.plot(t, S, 'b', alpha=0.7, linewidth=2, label='Susceptible')\n", + " #ax.plot(t, E, 'y', alpha=0.7, linewidth=2, label='Exposed')\n", + " ax.plot(t, I, 'r', alpha=0.7, linewidth=2, label='Infected')\n", + " ax.plot(t, R, 'g', alpha=0.7, linewidth=2, label='Recovered')\n", + " if D is not None:\n", + " ax.plot(t, D, 'k', alpha=0.7, linewidth=2, label='Dead')\n", + " ax.plot(t, S+E+I+R+D, 'c--', alpha=0.7, linewidth=2, label='Total')\n", + " else:\n", + " ax.plot(t, S+E+I+R, 'c--', alpha=0.7, linewidth=2, label='Total')\n", + "\n", + " ax.set_xlabel('Time (days)')\n", + "\n", + " ax.yaxis.set_tick_params(length=0)\n", + " ax.xaxis.set_tick_params(length=0)\n", + " ax.grid(b=True, which='major', c='w', lw=2, ls='-')\n", + " legend = ax.legend(borderpad=2.0)\n", + " legend.get_frame().set_alpha(0.5)\n", + " for spine in ('top', 'right', 'bottom', 'left'):\n", + " ax.spines[spine].set_visible(False)\n", + " if L is not None:\n", + " plt.title(\"Lockdown after {} days\".format(L))\n", + " plt.show();\n", + "\n", + " if R0 is not None or CFR is not None:\n", + " f = plt.figure(figsize=(12,4))\n", + "\n", + " if R0 is not None:\n", + " # sp1\n", + " ax1 = f.add_subplot(121)\n", + " ax1.plot(t, R0, 'b--', alpha=0.7, linewidth=2, label='R_0')\n", + "\n", + " ax1.set_xlabel('Time (days)')\n", + " ax1.title.set_text('R_0 over time')\n", + " # ax.set_ylabel('Number (1000s)')\n", + " # ax.set_ylim(0,1.2)\n", + " ax1.yaxis.set_tick_params(length=0)\n", + " ax1.xaxis.set_tick_params(length=0)\n", + " ax1.grid(b=True, which='major', c='w', lw=2, ls='-')\n", + " legend = ax1.legend()\n", + " legend.get_frame().set_alpha(0.5)\n", + " for spine in ('top', 'right', 'bottom', 'left'):\n", + " ax.spines[spine].set_visible(False)\n", + "\n", + " if Alpha is not None:\n", + " # sp2\n", + " ax2 = f.add_subplot(122)\n", + " ax2.plot(t, Alpha, 'r--', alpha=0.7, linewidth=2, label='alpha')\n", + "\n", + " ax2.set_xlabel('Time (days)')\n", + " ax2.title.set_text('fatality rate over time')\n", + " # ax.set_ylabel('Number (1000s)')\n", + " # ax.set_ylim(0,1.2)\n", + " ax2.yaxis.set_tick_params(length=0)\n", + " ax2.xaxis.set_tick_params(length=0)\n", + " ax2.grid(b=True, which='major', c='w', lw=2, ls='-')\n", + " legend = ax2.legend()\n", + " legend.get_frame().set_alpha(0.5)\n", + " for spine in ('top', 'right', 'bottom', 'left'):\n", + " ax.spines[spine].set_visible(False)\n", + "\n", + " plt.show();" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Create SEIR PDEs\n", + "\n", + "def deriv(y, t, N, beta, gamma, delta, alpha, rho):\n", + " S, E, I, R, D = y\n", + " dSdt = -beta * S * I / N\n", + " dEdt = beta * S * I / N - delta * E\n", + " dIdt = delta * E - (1 - alpha) * gamma * I - alpha * rho * I\n", + " dRdt = (1 - alpha) * gamma * I\n", + " dDdt = alpha * rho * I\n", + " return dSdt, dEdt, dIdt, dRdt, dDdt" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Load infections, deaths and recover data into dataframes\n", + "\n", + "import pandas as pd\n", + "\n", + "africa_infections_url = \"https://raw.githubusercontent.com/dsfsi/covid19africa/master/data/time_series/africa_daily_time_series_cases.csv\"\n", + "africa_deaths_url = \"https://raw.githubusercontent.com/dsfsi/covid19africa/master/data/time_series/africa_daily_time_series_deaths.csv\"\n", + "africa_recoveries_url = \"https://raw.githubusercontent.com/dsfsi/covid19africa/master/data/time_series/africa_daily_time_series_recovered.csv\"\n", + "\n", + "africa_infections_df = pd.read_csv(africa_infections_url).transpose()\n", + "africa_infections_df.columns = africa_infections_df.iloc[0]\n", + "africa_infections_df = africa_infections_df.drop(\"Country/Region\")\n", + "africa_infections_df = africa_infections_df.drop([\"Lat\",\"Long\"])\n", + "\n", + "africa_deaths_df = pd.read_csv(africa_deaths_url).transpose()\n", + "africa_deaths_df.columns = africa_deaths_df.iloc[0]\n", + "africa_deaths_df = africa_deaths_df.drop(\"Country/Region\")\n", + "africa_deaths_df = africa_deaths_df.drop([\"Lat\",\"Long\"])\n", + "\n", + "africa_recoveries_df = pd.read_csv(africa_recoveries_url).transpose()\n", + "africa_recoveries_df.columns = africa_recoveries_df.iloc[0]\n", + "africa_recoveries_df = africa_recoveries_df.drop(\"Country/Region\")\n", + "africa_recoveries_df = africa_recoveries_df.drop([\"Lat\",\"Long\"])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Country/RegionAlgeriaAngolaBeninBurkina FasoCabo VerdeCameroonCARChadCongoDRC...UgandaZambiaZimbabweBotswanaBurundiSierra LeoneMalawiSouth SudanWestern SaharaSao Tome and Principe
2020-01-220000000000...0000000000
2020-01-230000000000...0000000000
2020-01-240000000000...0000000000
2020-01-250000000000...0000000000
2020-01-260000000000...0000000000
\n", + "

5 rows × 53 columns

\n", + "
" + ], + "text/plain": [ + "Country/Region Algeria Angola Benin Burkina Faso Cabo Verde Cameroon CAR Chad \\\n", + "2020-01-22 0 0 0 0 0 0 0 0 \n", + "2020-01-23 0 0 0 0 0 0 0 0 \n", + "2020-01-24 0 0 0 0 0 0 0 0 \n", + "2020-01-25 0 0 0 0 0 0 0 0 \n", + "2020-01-26 0 0 0 0 0 0 0 0 \n", + "\n", + "Country/Region Congo DRC ... Uganda Zambia Zimbabwe Botswana Burundi \\\n", + "2020-01-22 0 0 ... 0 0 0 0 0 \n", + "2020-01-23 0 0 ... 0 0 0 0 0 \n", + "2020-01-24 0 0 ... 0 0 0 0 0 \n", + "2020-01-25 0 0 ... 0 0 0 0 0 \n", + "2020-01-26 0 0 ... 0 0 0 0 0 \n", + "\n", + "Country/Region Sierra Leone Malawi South Sudan Western Sahara \\\n", + "2020-01-22 0 0 0 0 \n", + "2020-01-23 0 0 0 0 \n", + "2020-01-24 0 0 0 0 \n", + "2020-01-25 0 0 0 0 \n", + "2020-01-26 0 0 0 0 \n", + "\n", + "Country/Region Sao Tome and Principe \n", + "2020-01-22 0 \n", + "2020-01-23 0 \n", + "2020-01-24 0 \n", + "2020-01-25 0 \n", + "2020-01-26 0 \n", + "\n", + "[5 rows x 53 columns]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Change date format\n", + "\n", + "africa_infections_df.index = pd.to_datetime(africa_infections_df.index)\n", + "africa_deaths_df.index = pd.to_datetime(africa_deaths_df.index)\n", + "africa_recoveries_df.index = pd.to_datetime(africa_recoveries_df.index)\n", + "\n", + "africa_recoveries_df.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Load population data\n", + "\n", + "pop_url = \"https://raw.githubusercontent.com/datasets/population/master/data/population.csv\"\n", + "pop_df = pd.read_csv(pop_url)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# install widgets\n", + "\n", + "#!pip install ipywidgets\n", + "#!jupyter nbextension enable --py widgetsnbextension" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d4341689d82c4d59a192b9a54e47c780", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(Dropdown(description='country', options=('Algeria', 'Angola', 'Benin', 'Burkina Faso', '…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "import ipywidgets as widgets\n", + "from ipywidgets import interact, interact_manual\n", + "import datetime\n", + "\n", + "@interact \n", + "def plot_SEIR(country = africa_infections_df.columns,\n", + " date = widgets.DatePicker(value=pd.to_datetime(africa_infections_df.index[-1])),\n", + " D = widgets.FloatText(value= 14), # how many days infection lasts for\n", + " incubation_period = widgets.FloatText(value= 1), # number of days an exposed individual becomes infectious\n", + " alpha = widgets.FloatText(value=0.002), # death rate\n", + " no_days_to_death = widgets.FloatText(value=14) \n", + " ):\n", + " \n", + " date = date.strftime(\"%Y-%m-%d\")\n", + " infections = africa_infections_df.loc[date, country]\n", + " no_deaths = africa_deaths_df.loc[date, country]\n", + " no_recoveries = africa_recoveries_df.loc[date, country]\n", + " R_0 = infections/no_recoveries\n", + "\n", + " N = pop_df.loc[(pop_df[\"Country Name\"] == country) & (pop_df[\"Year\"] == 2018)][\"Value\"].values\n", + " \n", + " gamma = 1.0 / D\n", + " delta = 1.0 / incubation_period # incubation period of five days\n", + " beta = R_0 * gamma # R_0 = beta / gamma, so beta = R_0 * gamma\n", + " rho = 1/no_days_to_death # days from infection until death \n", + " S0, E0, I0, R0, D0 = N-15515, 264, 15515, 7006, 264 # initial conditions: one exposed\n", + " \n", + " t = np.linspace(0, 600, 150) # Grid of time points (in days)\n", + " y0 = S0, E0, I0, R0, D0 # Initial conditions v\n", + " \n", + " ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta, alpha, rho))\n", + " S, E, I, R, D = ret.T\n", + " plotseird(t, S, E, I, R, D)\n", + " \n", + " return" + ] + } + ], + "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.8.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}