comprendre_SIR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SEIQR pour essayer de comprendre la dynamique de l'épidémie covid-19 (en France)\n",
    "\n",
    "_Disclamer: Je ne suis pas épidémiologiste._\n",
    "\n",
    "## Modèle\n",
    "Le modèle suppose 5 compartiments:\n",
    " - $S$: Compartiment des Succeptibles.\n",
    " - $E$: Compartiment des exposés (les personnes en incubation).\n",
    " - $I$: Compartiment des inféctés contagieux qui contaminent.\n",
    " - $Q$: Compartiemnt des infectés en quarantaine.\n",
    " - $R$: Compartiment des guéris et des morts (qui sont confondus dans ce\n",
    "   modèle).\n",
    "\n",
    "## Dynamique\n",
    "La dynamique du modèle est régie pas les équations suivantes:\n",
    " - $\\frac{dS(t)}{dt} = -\\beta\\frac{S(t)I(t)}{N}$\n",
    " - $\\frac{dE(t)}{dt} = \\beta\\frac{S(t)I(t)}{N} - aE(t)$\n",
    " - $\\frac{dI(t)}{dt} = aE(t)-\\gamma I(t)$\n",
    " - $\\frac{dQ(t)}{dt} = \\gamma I(t) - \\delta Q(t)$\n",
    " - $\\frac{dR(t)}{dt} = \\delta Q(t)$\n",
    "\n",
    "Sur les paramètres, on peut dire que:\n",
    " - $R_0=\\frac{\\beta}{\\gamma}$, et ça c'est interessant, car on connait $R_0$\n",
    "   (enfin on est capable d'aller dans la littérature chercher des valeurs\n",
    "   plausibles).\n",
    " - $a^{-1}$ est le temps moyen d'incubation.\n",
    " - $\\gamma^{-1}$ est le temps moyen pendant lequel un individu infectieux\n",
    "   contamine.\n",
    " - $\\delta^{-1}$ est le temps moyen qu'un indivudu contagieux qui ne contamine\n",
    "   plus (car isolé) reste en traitement avant mort ou guérison.\n",
    "\n",
    "## Initialisation\n",
    " - $I(0)$: nombre d'infectés non guéris et non morts (oui de considérer que tous le monde est contagieux dès le départ,\n",
    "   cela induit peut être le fait de voir le pic 1 ou 2 jours trop tôt).\n",
    " - $E(0)$ de telle manière à avoir $aE(0)$ qui correspond au flux (en utilisant\n",
    "   les données de la veille).\n",
    " - $Q(0)=0$, puisque j'ai mis tous les infectés dans $I$.\n",
    " - $R(0)=$, guéris + morts.\n",
    " - $S(0) = N-E(0)-I(0)-Q(0)-R(0)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Le code\n",
    "\n",
    "Cette partie ne vous interesse pas forcément. Si vous voulez juste utiliser, executez les cellules une à une, et rendez-vous à la partie utilisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ce dont on a besoin\n",
    "import collections\n",
    "from typing import List\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as tic\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Les types dont on a besoin\n",
    "State = collections.namedtuple(\"State\", (\"t\", \"S\", \"E\", \"I\", \"Q\", \"R\"))\n",
    "\n",
    "SimulationParam = collections.namedtuple(\n",
    "    \"SimulationParam\",\n",
    "    (\"incubation_period\", \"infection_period\", \"quarantaine_period\", \"R0\"),\n",
    ")\n",
    "\n",
    "InitializationParam = collections.namedtuple(\n",
    "    \"InitializationParam\",\n",
    "    ( \n",
    "        \"infected_today\",\n",
    "        \"infected_yesterday\",\n",
    "        \"recovered_today\",\n",
    "        \"recovered_yesterday\",\n",
    "        \"population_size\",\n",
    "    ),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La fonction principale de simulation\n",
    "def simu(\n",
    "    sim_par: SimulationParam, init_par: InitializationParam, timelimit: int = 365 * 24\n",
    ") -> List[State]:\n",
    "    \"\"\"Perform a simulation\n",
    "\n",
    "    Parameters:\n",
    "      - sim_par: A simulation parameter, all periods should be provided in days\n",
    "      - init_par: A initialization parameter, death are counted with recovered\n",
    "      - timeelimit: in hour, default (one year)\n",
    "\n",
    "    Returns:\n",
    "      - List of state, time in step is in hour.\n",
    "    \"\"\"\n",
    "\n",
    "    a = 1 / (sim_par.incubation_period * 24)\n",
    "    γ = 1 / (sim_par.infection_period * 24)\n",
    "    δ = 1 / (sim_par.quarantaine_period * 24) if sim_par.quarantaine_period > 0 else 1 \n",
    "    β = sim_par.R0 * γ \n",
    "\n",
    "    s_E = ( \n",
    "        init_par.infected_today\n",
    "        - init_par.infected_yesterday\n",
    "        + init_par.recovered_today\n",
    "        - init_par.recovered_yesterday\n",
    "    ) / (24 * a)\n",
    "    s_I = init_par.infected_today\n",
    "    s_R = init_par.recovered_today\n",
    "    s_Q = 0 \n",
    "    state = State(\n",
    "        0, init_par.population_size - s_E - s_I - s_R - s_Q, s_E, s_I, s_R, s_Q\n",
    "    ) \n",
    "\n",
    "    states = [state]\n",
    "    for t in range(1, 1 + timelimit):\n",
    "        succ_to_exposed = β * state.S * state.I / init_par.population_size\n",
    "        exposed_to_infected = a * state.E\n",
    "        infected_to_quarantaine = γ * state.I\n",
    "        quarantaine_to_recovered = δ * state.Q\n",
    "\n",
    "        state = State(\n",
    "            t=t,\n",
    "            S=state.S - succ_to_exposed,\n",
    "            E=state.E + succ_to_exposed - exposed_to_infected,\n",
    "            I=state.I + exposed_to_infected - infected_to_quarantaine,\n",
    "            Q=state.Q + infected_to_quarantaine - quarantaine_to_recovered,\n",
    "            R=state.R + quarantaine_to_recovered,\n",
    "        )\n",
    "\n",
    "        states.append(state)\n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La fonction de dessin\n",
    "def do_plots(\n",
    "    today: datetime.datetime,\n",
    "    sim_par: SimulationParam,\n",
    "    init_par: InitializationParam,\n",
    "    ax=plt,\n",
    "    label=None,\n",
    "    to_plot=lambda x: x.I + x.Q,\n",
    "    horizon=365,\n",
    "):\n",
    "    \"\"\"Plot a simulation\n",
    "\n",
    "    Parameters:\n",
    "      - today: datetime object describing the time t=0.\n",
    "      - sim_par: A simulation parameter, all periods should be provided in days\n",
    "      - init_par: A initialization parameter, death are counted with recovered\n",
    "      - ax: ax to plot (by default plt)\n",
    "      - label\n",
    "      - to_plot: a lambda function extracting fields to plot from a state\n",
    "      - horizon: in days\n",
    "    \"\"\"\n",
    "\n",
    "    sim = simu(sim_par, init_par, timelimit=horizon*24)\n",
    "    ax.plot(\n",
    "        [today + datetime.timedelta(x.t / 24) for x in sim],\n",
    "        [to_plot(x) for x in sim],\n",
    "        label=label,\n",
    "    )\n",
    "    if not hasattr(ax, 'yaxis'):\n",
    "        ax = ax.gca()\n",
    "    ax.yaxis.set_major_formatter(tic.EngFormatter(sep=\"\"))\n",
    "    ax.grid(True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utilisation\n",
    "\n",
    "Il convient dans un premier temps de définir ce qui va nous servir pour initialiser le modèle. J'utilise pour cela les données du 3 mars en France."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "today = datetime.datetime(2020,3,3)\n",
    "ip = InitializationParam(\n",
    "    infected_today=212, \n",
    "    infected_yesterday=191, \n",
    "    recovered_today=16, \n",
    "    recovered_yesterday=15, \n",
    "    population_size=67e6\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On fait tout d'abort un premier petit graphe avec une condition de simulation simple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD5CAYAAADLL+UrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAAEx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAucG9zdDE1ODkrZ2RmYTgzNjc3ZiwgaHR0cDovL21hdHBsb3RsaWIub3JnL4oVeZ0AACAASURBVHic7d17fFxnfefxzyON7ndbsqTIjpXEueKQi00gBFhlQyCw7abQliTbpYRN69KWV9rXQruQLm12u93SFlJoWbqkJIS2bNIutCVAQpIGi1zoEuzgxCbGl/gu27JsWaO7NJr57R/njD2WdZmRZuacmfm+Xy+9NHPmnDnPT7K/8+g5zznHmRkiIlIayoJugIiI5I9CX0SkhCj0RURKiEJfRKSEKPRFREqIQl9EpIREFlvBObcG+BugA0gAD5rZ551z9wO/Cgz4q95nZk8453qAzcCvmNlD/ntcB7wM/I6ZfWah/bW2tlp3d3dGRYyNjVFXV5fRNmGlWsKrWOopljqgOGrJRg1bt249aWZtaa1sZgt+AZ3A9f7jBmA3cBVwP/DxOdbvAV4Fnk5Z9ifAtrnWn/21YcMGy9TmzZsz3iasVEt4FUs9xVKHWXHUko0agC22SLYmvxYd3jGzY2b2sv94BNgJdC2y2SGg2jnX7pxzwG3Ak2l9ComISM44y+CMXOdcN/AcsB74z8DdwDCwBfiYmZ32h3c+DjyNNxz0Y+BXgIPAqM0xvOOc2wRsAmhvb9/w2GOPZVTE6Ogo9fX1GW0TVqolvIqlnmKpA4qjlmzUcPPNN281s41prZzunwRAPbAVeL//vB0oxzsY/EfAw3Z2eOfbeMcAvg98AXgX8wwHzf7S8M7moJuQNcVUi1nx1FMsdZgVRy2hG94BcM5VAN8AvmZm/+h/WPSbWdzMEsBfAzfM+jA5DsSAW4Fn0/oEEhGRnEpn9o4DHgJ2mtkDKcs7zeyY//R9wI45Nv99YJWZxb23ERGRIC0a+sBNwAeB7c65bf6y+4C7nHPXAgYcAH5t9oZm9oMstVNERLJg0dA3sxeAubrpT8yzfi/QO8fy+zNrmoiIZJvOyA25RML4xtYj9O46EXRTRKQIpDO8IwF6+MX9/I/v7ATg7+55M2+7tDXgFolIIVNPP8TiCeNLz+1jw9oWuppr+MLmPUE3SUQKnEI/xLYdHmJgZIpfvnEtd75pDf9v3yAnhieDbpaIFDCFfog9t3uAMgc9l63inVe1A9C7a2CRrURE5qfQD7HtfVEuXdVAU20FV3Q00NZQxQ9ePxl0s0SkgCn0Q8rMePVIlPVdTQA457j+wma2HR4KuGUiUsgU+iHVPzzFydEpru5qPLPsugtbOHBqnMGx6QBbJiKFTKEfUrv6RwC4svNs6F+7phmAV46oty8iS6PQD6kDJ8cAuKjt7B11ruzwPgB2Hx8JpE0iUvgU+iF14NQYdZXltNVXnVnWVFtBe2PVmb8CREQypdAPqQMnx1i7so7ZVye9rL2B3Qp9EVkihX5IHTw1Tndr7XnLL29vYE//KPFE+nc8ExFJUuiH0Ew8weHT43SvrDvvtcvaG5iaSXBocDyAlolIoVPoh9CJkSlicWN1y/k9/UtWeR8EyQO9IiKZUOiH0HH/+jqdTdXnvbbW7/0fOKXQF5HMKfRD6HjUC/32xvNDf2VdJfVVEQ6e0vCOiGROoR9CydCfq6fvnGPtylr19EVkSRT6IXR8eJLKSBnNtRVzvt69sk49fRFZEoV+CB2PTtLZVH3eHP2ktStrOTw4zkw8keeWiUihU+iH0PHo5Jzj+UndK+uYSRhHh3RDFRHJjEI/hI4PT845np+0dqU3lVPj+iKSKYV+yJgZx4cn6Vigp796hRf6R05P5KtZIlIkFPohMxaD6ZkEqxYI/faGKsrLHH1DOpgrIplR6IfM8LR3TZ3W+sp514mUl9HRWE2fevoikiGFfsiMnAn9qgXX62qpoW9IoS8imVHoh0yyp79ygZ4+wOrmGvX0RSRjCv2QORP6dYv39I8PTxLTXH0RyYBCP2SGp7zQb5nnbNykruYaEnb2kg0iIulQ6IfMyLTRUltBpHzhX01XSw2AxvVFJCMK/ZAZnjZWLnIQF7yePqBxfRHJiEI/ZEamjZV1Cx/EBbigWT19EcncoqHvnFvjnNvsnNvpnPuJc+63/OUrnHPPOOf2+N9b/OV3O+fMOXdLynu8z1/2C7krpTgMT9ui0zUBqivKaa2v4qhCX0QykE5Pfwb4mJldCbwF+E3n3FXAJ4BnzexS4Fn/edJ24K6U53cCr2SnycVteNpYkUZPHzRXX0Qyt2jom9kxM3vZfzwC7AS6gNuBr/qrfRX4uZTNngducM5VOOfqgXXAtmw2vBjF4gnGYovP0U9a3azQF5HMRDJZ2TnXDVwH/BBoN7Nj4H0wOOdWpaxqwL8A7waagMeBixZ4303AJoD29nZ6e3szaRajo6MZbxNGQ5PenPvBowfp7T266Pqx4Sn6BmfYvHnzvNfeD1Kx/F6SiqWeYqkDiqOWfNeQduj7PfZvAL9tZsNphMxjwL14of8x4L75VjSzB4EHATZu3Gg9PT3pNguA3t5eMt0mjHYeG4be57nxuvX0XN256Pp7yvbx9MGdbLjxbTRWLzyvPwjF8ntJKpZ6iqUOKI5a8l1DWrN3nHMVeIH/NTP7R39xv3Ou03+9EziRuo2ZvQSsB1rNbHf2mly8To1OA6Q9pt/hX3NfJ2iJSLrSmb3jgIeAnWb2QMpLjwMf8h9/CPjmHJt/kgV6+HKuoQkv9Jtr0wv95I1Wjin0RSRN6Qzv3AR8ENjunEsejL0P+DTwD865e4BDwC/O3tDMnsxWQ0tBdCIGMO8N0Wc729PXwVwRSc+ioW9mLwDzDeDfMnuBmT0CPDLH8rsza1rpGRr3Qr+pJr3QX9VQjXPq6YtI+nRGbogMT8SIlHknXqWjMlLGyroqjemLSNoU+iESnYhRX5HZ1MvOpmqODyv0RSQ9Cv0QGRqPkeZw/hkdTdXq6YtI2hT6IRKdiFEXybynrzF9EUmXQj9EohMx6jIc3uloqiY6EWN8eiZHrRKRYqLQD5GlhH6nTtASkQwo9EMkOpH5mH57o0JfRNKn0A+JWDzB6NTMEmbveDdT0bi+iKRDoR8Sw/7ZuLWZjukne/qatikiaVDoh0TyEgyZjunXVJbTXFuh4R0RSYtCPyTOhn7m23Y0atqmiKRHoR8SQ8nQz3CePvgnaA3romsisjiFfkgsdUwf/EsxqKcvImlQ6IdEcngn09k7AB2NNZwcnWZqJp7tZolIkVHoh0R0PNnTz3zb5AlaJ4anstkkESlCCv2QGJqIUVtZTqRsaWP6oGmbIrI4hX5IRCdiNKd585TZdNtEEUmXQj8kohMxGpcY+rptooikS6EfEtHxWNq3SZytobqC+qqIevoisiiFfkhEJ2Jp3xB9LrqZioikQ6EfEtGJpff0QTdTEZH0KPRDYmhielmh39Gonr6ILE6hHwJTM3EmY4nlhX5TNSdGJpmJJ7LYMhEpNgr9EEiejdtUW7nk9+hoqiZhMDCqE7REZH4K/RBIno273DF90Fx9EVmYQj8EzvT0lzWm791BS+P6IrIQhX4IJEN/qWfkgnr6IpIehX4IDGVheKe5toKqSJnOyhWRBSn0QyAbwzvOOc3VF5FFKfRDIBn6S732TpLOyhWRxSj0QyA6EaOhOkL5Ei6rnKqzqUY9fRFZkEI/BJZ7CYakjqZq+ocnSSQsC60SkWK0aOg75x52zp1wzu1IWXa/c67PObfN/3qvv7zHOWfOuXtS1r3OX/bx3JRQ+JZ7sbWkzqZqZhLGyTGdoCUic0unp/8IcNscy//czK71v55IWb4duCPl+Z3AK0tvYvEbGl/edXeSOhqT19XXEI+IzG3R0Dez54DBDN7zEFDtnGt3zjm8D4wnl9i+kpCt4Z3OJu8ELY3ri8h8IsvY9qPOuV8GtgAfM7PTKa99HfhF4MfAy8CC4w3OuU3AJoD29nZ6e3szasjo6GjG24TJQHScNVVT9Pb2LquW6JQ3lv/8lu1UDfw0iy1cmkL/vcxWLPUUSx1QHLXkvQYzW/QL6AZ2pDxvB8rx/lL4I+Bhf3kP8G2gA/g+8AXgXcD9wMfT2deGDRssU5s3b854m7BIJBJ26X1P2B8/sdPMlldLPJ6wdfd9xz795M4stW55Cvn3MpdiqadY6jArjlqyUQOwxdLIVzNb2uwdM+s3s7iZJYC/Bm6Y9fpxIAbcCjy7lH2UiolYnOn48i6rnFRW5mjXdfVFZAFLGt5xznWa2TH/6fuAHXOs9vvAKjOLe0P7Mpcz193JwuwdSN5BS5diEJG5LRr6zrlH8YZtWp1zR4A/AHqcc9cCBhwAfm32dmb2g6y2tEhl4xIMqdobq9nRF83Ke4lI8Vk09M3srjkWPzTPur1A7xzL78+wXSUjGxdbS9XZVM0zr/VjZugvLBGZTWfkBizbPf2OphqmZhJnPkxERFIp9AOW7dDXdfVFZCEK/YCduVVilg7kdvihf3xYB3NF5HwK/YBFJ2KUlzkaqpZzntxZ6umLyEIU+gGLTsRorI5k7aBrW30VZU7X3xGRuSn0AzaUpevuJEXKy1jVoDtoicjcFPoBi07EaKqtzOp76g5aIjIfhX7AsnWFzVQ6K1dE5qPQD1g0S9fST9Xh3yDduw6TiMhZCv2AeT397MzcSepsqmZ8Os7I1ExW31dECp9CP0CJhHm3SqzJ9pi+fzOVIY3ri8i5FPoBGp2eIWHZOxs3qavZC/2+ofGsvq+IFD6FfoCyfTZu0uoWP/RP62CuiJxLoR+gbF93J6mtvorK8jKODCn0ReRcCv0AnbmBSpZDv6zMcUFzNUfU0xeRWRT6AUpe/rg5yydnAXS11Gh4R0TOo9APUK6Gd8A7mNun4R0RmUWhH6ChiWkge/fHTdXVXMvAyBSTsXjW31tECpdCP0DR8RiVkTKqK8qz/t7JGTxH1dsXkRQK/QB5J2Zlv5cP3pg+oCEeETmHQj9AQ+OxnAztQMoJWjqYKyIpFPoBysUVNpM6m6opL3Oatiki51DoB8i7gUr2p2uCdzOVjsZqDe+IyDkU+gHKxWWVU3U1a66+iJxLoR+g6ETuxvTBm8Gjnr6IpFLoByQWTzA2Hc/Z7B3wZvAci04Qiydytg8RKSwK/YCcORs3hz39ruYaEobulysiZyj0A5K87k5Ox/Q1V19EZlHoByTqX4Ihl6G/uqUWgMODupmKiHgU+gE5c1nlHFxhM6mruYYyB4cU+iLiU+gH5MxllXPY06+MlHFBcw0HTyn0RcSj0A9IPsb0AbpX1nFQPX0R8S0a+s65h51zJ5xzO1KWrXDOPeOc2+N/b/GX3+2cM+fcLSnrvs9f9gu5KaEwJYd3GnMc+heurOXQqbGc7kNECkc6Pf1HgNtmLfsE8KyZXQo86z9P2g7clfL8TuCVZbSxKEUnYjRURygvczndz9oVtZwej535kBGR0rZo6JvZc8DgrMW3A1/1H38V+LmU154HbnDOVTjn6oF1wLYstLWo5Pps3KS1K+sAOKRxfREBnJktvpJz3cC3zWy9/3zIzJpTXj9tZi3OubuBjcA08D2gCS/0L/K3//o8778J2ATQ3t6+4bHHHsuoiNHRUerr6zPaJmgPbJ1keMq4/6015yzPdi2HRxJ86sUJfuOaKm7ojGTtfdNRiL+XhRRLPcVSBxRHLdmo4eabb95qZhvTWTdXKfAYcC9e6H8MuG+hlc3sQeBBgI0bN1pPT09GO+vt7SXTbYL2+ddeZHVThJ6eN5+zPNu1jE3N8KkXn6Kuo5uennVZe990FOLvZSHFUk+x1AHFUUu+a1jq7J1+51wngP/9ROqLZvYSsB5oNbPdy2ticcrlDVRS1VVFaK2v0vCOiABLD/3HgQ/5jz8EfHOOdT7JIj38UjY4Ns2KutydmJWqe2UtBwc1g0dE0puy+Sjwr8Dlzrkjzrl7gE8Dtzrn9gC3+s/PYWZPmtnmbDe4GMzEE0QnYnkL/QtX1uoELREB0hjTN7O75nnpltkLzOwRvCmes5ffnWG7itqQP30yX6G/dkUd//TjPiZjcaoryvOyTxEJJ52RG4DTY97F1lpyeN2dVGtX1mIGR06rty9S6hT6ARj0Qz9vPf2V3tU2D5xU6IuUOoV+AAbz3NO/uM2bA7zv5Ghe9ici4aXQD8DgeH57+k01FbQ1VLH3hEJfpNQp9AOQHNPPxzz9pEva6nh9QNM2RUqdQj8Ag2Mx6irL8zqT5pK2evaeGCWdy26ISPFS6Afg9Pg0LXka2km6pK2e6ETszPEEESlNCv0ADI5NszLfob/KO5irIR6R0qbQD0AwPX3vEsuvD+hgrkgpU+gHYHBsmhV5mq6ZdEFTDdUVZbyuGTwiJU2hH4DBsfz39MvKHBe31qunL1LiFPp5NhmLMz4dz9sc/VSXrKrXmL5IiVPo59np8fyejZvqkrY6Dp8eZzIWz/u+RSQcFPp5dmo0v2fjprp0VQNm6MxckRKm0M+zgdEpANoaqvK+7ys6GwD46fGRvO9bRMJBoZ9nAyNe6K8KIPS7V9ZRFSnjp8eG875vEQkHhX6enfR7+q31+Q/98jLH5R0N7OpXT1+kVCn082xgZIr6qgg1lcHcwery9gZ2HlPoi5QqhX6eDYxM0Vqf/4O4SVd0NnJydOrMMJOIlBaFfp6dHJ0K5CBu0pUd3sHcXTqYK1KSFPp5NjASbOhf3pGcwaODuSKlSKGfZydHpwM5iJu0sr6KVQ1VGtcXKVEK/TyamokTnYjRFmDog9fb36lpmyIlSaGfRyf9s3FbAxzeAVjf1cTu/hFdjkGkBCn08+ikP2Mm6J7+NaubmEmYevsiJUihn0fJaZJBHsgFeOPqZgC290UDbYeI5J9CP4+CvO5Oqs6malrrK3nlsEJfpNQo9PPoWHSSMhd86DvneOPqZl49MhRoO0Qk/xT6eXQ8OkFbQxUV5cH/2N+4uom9A6OMTs0E3RQRyaPg06eEHItO0tFUE3QzALhmdTNm8BON64uUFIV+Hh2LTtLZWB10MwC4enUTAK9oiEekpCj08+h4dJKOpnCEfmt9FReuqGXLgdNBN0VE8mhZoe+cO+Cc2+6c2+ac2+Ive8Q5N+6ca0hZ7/POOXPOtS63wYVqZDLG6NQMFzSHI/QB3tS9gi0HT2NmQTdFRPIkGz39m83sWjPbmLJsL3A7gHOuDLgZ6MvCvgrW8egkQGjG9AFuuKiFwbFpXh/QPXNFSkWuhnceBe7wH/cALwIlPU3kmB/6nSEZ3gGvpw/w0n4N8YiUisgytzfgaeecAV8yswf95XuA251zLcBdwN8B75nvTZxzm4BNAO3t7fT29mbUiNHR0Yy3ybfnjsQAOPDaNsYOzP9Zm89azIzGSvjWD1/jgol9WX//Qvi9ZKJY6imWOqA4asl7DWa25C/gAv/7KuAV4B3AI8AvAL8L/DrwKt5fFAeA1sXec8OGDZapzZs3Z7xNvj3w9C7r/sS3bSoWX3C9fNfykb/dYm/942dz8t6F8HvJRLHUUyx1mBVHLdmoAdhiaeb2soZ3zOyo//0E8E/ADSkvPwb8IfCMmSWWs59icPj0OJ2N1VRGwjVh6k3dK+gbmqBvaCLopohIHiw5gZxzdckZOs65OuBdwI7k62Z2CPg94IvLbWQxODw4zpoVtUE34zw3rfMmVL2wZyDglohIPiyn29kOvOCcewV4CfiOmX03dQUz+5KZvb6cBhaLQyEN/cva62lvrOK5PSeDboqI5MGSD+Sa2T7gmjmW3z3P+t1L3Vehm4zF6R+e4sIQhr5zjrdf2sYzr/UTTxjlZS7oJolIDoVrgLlIHTk9DhDK0Ad4x2VtRCdiuuqmSAlQ6OfB4UHvIGkYh3cA3rauFefgeQ3xiBQ9hX4eHBoMd09/RV0lV3c10bvrRNBNEZEcU+jnwcFT49RUlNNaXxl0U+b1zivbefnQECeGJ4NuiojkkEI/D/YOjHLJqjqcC+9B0ves7wDgqdf6A26JiOSSQj8P9vaPsK6tPuhmLGjdqnoubqvjqR3Hg26KiOSQQj/HxqZmOBqdZN2qcIe+c47b3tDBv+47xdD4dNDNEZEcUejnWPKyxetWNSyyZvBuW99BPGE89RP19kWKlUI/x/b0J0M/3D19gKu7mri4tY5vvFzStz4QKWoK/RzbOzBKRblj7cpwTtdM5Zzj5zes5qX9gxw8NRZ0c0QkBxT6Ofba0WEuaaunorwwftTvv74L51BvX6RIFUYSFSgzY0dflPVdTUE3JW2dTTW8bV0r39h6hHhC984VKTYK/Rw6PjzJqbFpri6g0Af4DzdcSN/QBM9ozr5I0VHo59D2I1GAgurpA9x6VTtdzTU8/OL+oJsiIlmm0M+hHX1Ryhxc1dkYdFMyEikv4+63dvPS/kF29EWDbo6IZJFCP4dePjTEZe0N1FSWB92UjN1xwxrqKsv5q+/rHjgixUShnyOxeIKtB0/zlotXBt2UJWmsruDDN13Ed149xs5jw0E3R0SyRKGfI9v7okzE4rz5ohVBN2XJfvXtF9NQHeGzT+8OuikikiUK/Rz54b5BAN5UwKHfVFvBprdfzL/s7Oel/YNBN0dEskChnyPP7xngsvZ6Wuurgm7Kstzz9ovoaq7hU/+8g1g8EXRzRGSZFPo5EB2P8cP9g9xyZXvQTVm22soIf/CzV7Grf4SvaAqnSMFT6OdA7+4TxBPGO4sg9MGbt//OK9v5zNO7+elxHdQVKWQK/Rx4cvtxWusruXZNc9BNyQrnHH/8/qtprK7g3kd/zGQsHnSTRGSJFPpZdmp0imd/2s/t13ZRXhbe2yNmqq2higc+cA27+0f53a+/ipmuyyNSiBT6WfbP244Sixsf2Lgm6KZk3Tsua+N33n05j79ylD9/RtM4RQpRJOgGFJOZeIJHfrCf6y5s5vKO8N8payl+o+cSDp0a5y++t5eqinJ+8+Z1QTdJRDKg0M+ix185yuHBCX7/Z94QdFNyxjnHH71vPVMzcf7sqV2MTc3w8XddTlkRDWWJFDOFfpaMTc3wmad2cVVnI7dcsSro5uRUpLyMz37gWmoqy/li7+vs7h/lgTuuobG6IuimicgiNKafJQ88s5uj0Un+++1vKIleb3mZ43++72ru/9mr2LzrBO/53PP07joRdLNEZBEK/Sx4cvsxHnphPx98y1o2dhfuZRcy5Zzj7psu4v9+5EZqKsu5+ys/4iN/u5U9/SNBN01E5qHhnWV6dmc/v/XYNq5d08x//Zkrg25OIK6/sIXv3Ps2vvT9fTz43D6efu04735DB1fXxPk3ZjhX/H/5iBQKhf4STcbifLH3db7wvT284YImHvnwm6iKFN5187OlKlLOvbdcyn98y1oefG4fj/3oEE+Ox3j09c2866oO3nVVO9dd2EJlRH9cigQp56HvnLsN+DxQDnzZzD6d633m0omRSR7fdpSHX9jP0egk77++iz+8fT11Vfr8BFhRV8kn3nMFv/3OS/nsP3yPvdP1/O2/HuShF/ZTFSnjmjXNbFjbwhUdDVzSVs8lbfUFeZMZkUKV06RyzpUD/wu4FTgC/Mg597iZvZbL/S6VmTE1k2BsaobRqRlGJmfoG5rgyOkJ9g2MsvXgaXb1j2AG11/YzJ/94jXctK416GaHUnVFOTd1VfB7PTcwOjXDC3sGeGn/abYeHOSvn9vHTMI7o9c5aKuvor2x2v+qorW+iobqCI01FTRWR2isrqC+OkJVpJzKSBlVkTIqk1/l3nMNIYmkJ9fd0xuAvWa2D8A59xhwO5D10P93f/E8E7E4ZpAwO+/72cdeuCfMMCCROLvOdDxBLD735QUaqiNcu6aZ917dyXuv7mDdquI8+SoX6qsi3La+k9vWdwIwNRPn4Klx9p4YZU//KEeHJjg+PMmR0+NsOTjI0Hgs431EyhxlzuGcN7No9uMyh//df1zmvT6b4/yFyfUmJiao+dHmOdZgzg+d85bMub/F3yfbxsfGqX35+znfTz4UQy3JGlbUVvIPH7kx5/vLdeh3AYdTnh8B3jx7JefcJmATQHt7O729vRntZHR0lEYi1JVDcrakc1CG9x/b+c/dmeWA/7wsuRxHpCxCdQRqIu7M9xXVjraaMuoqwLkJoI8jr/VxJEd/q4yOjmZcf1gtVkstcE0ErmkFzvzBVEnCKpiYgfGYMTFjjM/AxIwxk4BYAmIJ/3EcZhJGLAFxw/vwxvtQTz72PtzN/7DH+6D3P/xn87oB5y08Y6YsQXlkKq3aZ1+aKJ0rFeXrakYNNQkiZRN52ltuFUMtyRqqZybz8n8/16E/V7flvH/bZvYg8CDAxo0braenJ6Od9Pb28ui9mW0TVr29vWRaf1gVUy1QPPUUSx1QHLXku4ZcT6U4AqReeWw1cDTH+xQRkXnkOvR/BFzqnLvIOVcJ3Ak8nuN9iojIPHI6vGNmM865jwJP4U3ZfNjMfpLLfYqIyPxyPrnczJ4Ansj1fkREZHE6PVJEpIQo9EVESohCX0SkhCj0RURKiLPZpw4GzDk3ABzMcLNW4GQOmhME1RJexVJPsdQBxVFLNmpYa2Zt6awYutBfCufcFjPbGHQ7skG1hFex1FMsdUBx1JLvGjS8IyJSQhT6IiIlpFhC/8GgG5BFqiW8iqWeYqkDiqOWvNZQFGP6IiKSnmLp6YuISBoU+iIipcS7y1B+v/Cusb8Z2An8BPgtf/kK4Blgj/+9xV/+S8Cr/tcPgGtS3us2YBewF/jEAvv8kP++e4APpSz/LvCK347/DZQXcC29/vbb/K9VhVgL0JBSwza8OcyfK/B/Z3f47/sT4E8LoI7vAkPAt2ct/6i/rQGtAf9OHgZOADsW2eecNS+1lpDV8BBefr0KfB2oX7T9mf7SsvEFdALX+48bgN3AVcCfJgsCPgH8if/4rSk/wPcAP/QflwOvAxcDlX7xV82xvxXAPv97i/84+X6N/ncHfAO4s4Br6QU2FsPvZdZ6W4F3FGo9wErgENDmr/dV4Jaw1uGvewvws5wf+tcB3cABlhb6WanFf/4O4HoWCMyFal5qLSGroTFlvQdY4IP8zHqZ/tJy8QV8E7gV7HG1gQAAA1JJREFU75OsM+UHu2uOdVuAPv/xjcBTKa99EvjkHNvcBXwp5fmXgLtmrVMBfAu4o1BrYZmhH6ZaUpZdinefZVeo9QBvAv4lZfkHgS+GtY6U13uYFfoprx1gCaGfrVpSlnWzcGAuWvNyawlJDQ74K+C/LNbewMf0nXPdeJ+4PwTazewYgP991Ryb3AM86T+e68brXXNss+B6zrmn8P7EGsH7E2lJwlAL8BXn3Dbn3Kecc3PdozgtIakFvND8e/P/ZS9VwPXsBa5wznU75yLAz3HubUTDVkdeLLOWdOW05jDU4Jz7CnAcuAL4y8XeLOc3UVmIc64eb0jlt81seLGMcs7djPdDe1ty0RyrzRUOC65nZu92zlUDXwP+Ld54XEZCUssvmVmfc67Bb8sHgb9Jo/mz2xaGWpLuxKtjyYKux8xOO+d+Hfh7IIE3rntxms1PbVe+6si5LNSS9q7mWJaVmsNSg5l92DlXjhf4dwBfWejNAuvpO+cq8H5gXzOzf/QX9zvnOv3XO/F638n13wh8GbjdzE75i+e88bpz7s1+b3ebc+7fz7deanvMbBLv/r23F2otZtbnfx8B/g9wQ6HW4r/3NUDEzLZmWkfY6jGzb5nZm83sRrxhgD0hriOnslTLfO+9JqWWj5DG//1iqMHM4nidip9ftPHLHZNb4hiYw+uBfm7W8j/j3AMhf+o/vhDvT+S3zlo/gnew7CLOHuB4wxz7WwHsxxtPa/EfrwDqOTsGF/F/aB8t0Foi+OOSeMcnvg58pBBrSXn908B/K/R/Z/5rq/zvLXgzki4Lax0p6/eQ5TH9bNWSsl03C4+HL1pzprWEpQa/HetS2vQZ4DOLtn+p/6GW84X3543hTTNKTst7L94sh2fxekHPpvyH+TJwOmXdLSnv9V68o+evA7+3wD7/k/+D3wt82F/WDvyIs1Pp/hKvZ1mItdThzXJJ1vJ5Mp9+GopaUl7bB1xR6P/O/OWPAq/5X5nOEAuijueBAWACr6f5bn/5vf7zGbze5pcDrOVR4BgQ89t0zzz7nLPmpdYSlhrwRmpeBLYDO/CGpxsXa78uwyAiUkICn70jIiL5o9AXESkhCn0RkRKi0BcRKSEKfRGREqLQFxEpIQp9EZES8v8BfM6FQ0bOGSwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sp = SimulationParam(\n",
    "        incubation_period=5.4, \n",
    "        infection_period=2, \n",
    "        quarantaine_period=10, \n",
    "        R0=3)\n",
    "do_plots(today, sp, ip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Puis un graphe qui bouge en fonction de $R_0$ et de la période d'incubation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c7209a1c3ae4d9b9afcefe74bc1e3e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=3.5, description='R0', max=5.0, min=1.0), FloatSlider(value=5.4, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_computer_is_slow = False\n",
    "@interact(\n",
    "    R0=widgets.FloatSlider(3.5,min=1,max=5,continuous_update=not my_computer_is_slow),\n",
    "    incub=widgets.FloatSlider(5.4,min=2,max=14,continuous_update=not my_computer_is_slow),\n",
    "    horizon=widgets.IntSlider(250,min=10,max=366,continuous_update=not my_computer_is_slow)\n",
    ")\n",
    "def a_plot(R0, incub, horizon):\n",
    "    fig,(ax1,ax2) = plt.subplots(2,1, figsize=(12,8))\n",
    "    sp_ref = SimulationParam(\n",
    "        incubation_period=5.4, \n",
    "        infection_period=2, \n",
    "        quarantaine_period=10, \n",
    "        R0=1.5)\n",
    "    sp = SimulationParam(\n",
    "        incubation_period=incub, \n",
    "        infection_period=2, \n",
    "        quarantaine_period=10, \n",
    "        R0=R0)\n",
    "    do_plots(today, sp_ref, ip, ax=ax1, label=f'$R_0={sp_ref.R0:.2f}, incub={sp_ref.incubation_period:.2f}$', horizon=horizon)\n",
    "    do_plots(today, sp, ip, ax=ax1, label=f'$R_0={sp.R0:.2f}, incub={sp.incubation_period:.2f}$', horizon=horizon)\n",
    "    ax1.legend()\n",
    "    ax1.set_ylabel('I(t)+Q(t)')\n",
    "    do_plots(today, sp_ref, ip, ax=ax2, label=f'$R_0={sp_ref.R0:.2f}, incub={sp_ref.incubation_period:.2f}$', to_plot=lambda x: x.I+x.Q+x.R, horizon=horizon)\n",
    "    do_plots(today, sp, ip, ax=ax2, label=f'$R_0={sp.R0:.2f}, incub={sp.incubation_period:.2f}$', to_plot=lambda x: x.I+x.Q+x.R, horizon=horizon)\n",
    "    ax2.legend()\n",
    "    ax2.set_ylabel('I(t)+Q(t)+R(t)')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avec ce précédent graphique, on peut constater plein de choses, par exemple:\n",
    "\n",
    " - $R_0$ a un influence sur la hauteur du pic, donc du nombre de personne à traiter (dont une partie devra être hospitalisé).\n",
    " - la péride d'incubation a très peu d'influence sur la hauteur des pics et sur l'équilibre atteint, mais ça retarde la progression.\n",
    " \n",
    "Ma conclusion (qui n'engage que moi): On comprend aussi que si $R_0$ n'est pas faible le système de santé n'arrivera pas à traiter en même temps tous les patients. Or notre moyen d'action est sur $\\beta$ et comme $R_0=\\frac{\\beta}{\\gamma}$, nous pouvons tenir $R_0$ le plus faible possible en limitant la propagation (diminution contacts, lavage des mains, isolation des contagieux…)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}