model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling of the impact of social tracing apps: a contact network model\n",
    "\n",
    "23th April 2020\n",
    "\n",
    "by Alexandre Variengien, Simon Dupouy, Guilhem Niot, Timothé Picavet and Pierre Marrec.\n",
    "\n",
    "contact : [alexandre.variengien@ens-lyon.fr]\n",
    "\n",
    "We wish to thank the _Soyez Prevenus_ team for proposing and organising this modelling challenge !\n",
    "\n",
    "**An online interactive version of this document is available [here](https://mybinder.org/v2/gh/GuilhemN/covid-19-tracing-app-modelling/master?filepath=model.ipynb)**\n",
    "### Context of creation\n",
    "\n",
    "This model was created during the \"Hack COVID19\" Hackathon, which took place between the 10th and the 12th of April 2020 and was organized by HEC and l'École Polytechnique.\n",
    "It was created to model the impact of the adoption of a social tracing app, like STOP COVID in France, on the spread of the pandemic.\n",
    "Its primary goal is to show the influence of parameters such as the proportion of people using the app in the global population or the effect of individual behaviour after receiving an exposure notification from the app on the epidemic.\n",
    "\n",
    "\n",
    "### Motivation\n",
    "\n",
    "We are students in computer science at ENS de Lyon, so we are not experts in the field of epidemiology. This model doesn't aim to provide precise predictions for the future of the pandemic.\n",
    "However, it is easily understandable and still reasonnably effective to get an overall sense of what we could expect from a social tracing app.\n",
    "More specifically, we believe it gives a good intuition about the impact a social tracing app may have on the epidemic depending on the adoption rate of the app.\n",
    "You may try changing the parameters to convince yourself that targeted quarantining of suspected cases can lead to a great limitation of the spread - if it is well adopted by the population.\n",
    "So feel free to play with our model!\n",
    "\n",
    "### Bibliography\n",
    "\n",
    "For this model we used the contact network epidemiology approach, as it has been widely studied and it can account for the discrepency of transmission between individuals.\n",
    "It has been discussed that this approach is suited to model an epidemic as SARS-CoV-1 in [1].\n",
    "Moreover, the structure of the chosen graph has a great impact on the final extent of the outbreak, as discussed in [2].\n",
    "We chose two models for the graph:\n",
    "- a model using an exponential law for the degrees distribution - which seems suited for a disease such as SARS in urban areas [3]\n",
    "- a method from [3] for building the network by stages: first we build highly connected households, then the connections between them.\n",
    "\n",
    "As for the epidemiological parameters of the disease, we adapted parameters from [4] to suit our model.\n",
    "Also, other models for the impact of a social tracing app have been proposed such as in [4] and [5].\n",
    "\n",
    "### The model\n",
    "\n",
    "Each vertex of the graph is an individual that can be : Healthy, Asymptomatic, Presymptomatic, Symptomatic, Dead or Recovered.\n",
    "At each step we go trough all the edges and with the probability given by these edges, there is a contact between the individuals.\n",
    "There is a parameterized proportion of close contact which are likely to propagate the virus and of being detected by the app, whereas further contact are less likely to propagate the virus but are not detected by the app. The last type of contact models evironnemental contamination by surfaces for example.\n",
    "\n",
    "The model implements several policies based on screening when an individual gets symptoms:\n",
    "- Warning after symptoms:\n",
    "    * If this policy is applied, the app sends notifications to the persons the infected individual met in the past 14 days as soon as the symptoms are detected. To be more realistic, in this case, non infected individuals may also send notifications at random to model other diseases with similar symptoms.\n",
    "    * If not, the individual passes a screening test and notifications are only sent if the test is positive (the result is only available after a predefined number of days, which is approximately 5 days for now in France).\n",
    "\n",
    "- Quarantine after notification:\n",
    "    * If this policy is applied, a proportion of individuals receiving a notification goes in quarantine directly. They pass a screening test and stop the quarantine if it turns out to be negative.\n",
    "    * If not, individuals receiving a notification wait until they get the result of their test before going in quarantine.\n",
    "\n",
    "### Limits of the model\n",
    "\n",
    " * Some parameters like the proportion of asymptomatic persons and their infectiosity are not well known. While our work is based on their latest evaluation, they may be reevaluated quite differently in a few days or weeks.\n",
    "\n",
    " * We don't take into account the eventual false positives of nasal tests. We consider here that there are none.\n",
    "\n",
    " * The variability of infectiousness seems to be really high: an infected person that will develop symptoms is really contagious during the 2 days before the symptoms appear and is less contagious the 4 days after the symptoms appear, according to [6]. This could moderate our conclusions about the validity window in the \"results\" part as it gives us 2 to 3 free days during the incubation period.\n",
    "\n",
    " * The geographical heterogeneity of a national territory is not taken into account, this is a simulation just for one small isolated community (the simulations presented below in the result part are done with 4000 individuals)\n",
    "\n",
    " * The community is also homogeneous : differences in sex, age or prior conditions that would affect the infection are not taken into account.\n",
    "\n",
    " * In the initial state, a proportion of the population (5%) is contaminated at the same time. In reality, at the end of the lockdown, the infected will have been contaminated at different instant, smoothing the pressure on healthcare systems.\n",
    " \n",
    " * The distribution of the time an individual spends as asymptomatic or symptomatic follows a geometrical law. In reality it seems that they are more close to a normal or lognormal distribution. However, beacause symptomatic individuals are quarantined and asymptomatic ones are not the most infectious, this bias should not impact the dynamic of the epidemic too much.\n",
    "\n",
    "### Key takeaways\n",
    "\n",
    "\n",
    " * Contaminations from presymptomatic individuals are critical and social distancing can be effective to limit their effect. However, individuals must be warned early (since they do not yet display symptoms). The duration of the signaling chain must be less than the 5 days incubation period.\n",
    " * Even with quicker and more precise screening, the app seems unlikely to be effective because of the validity window of nasal tests.\n",
    " * The advice of quarantining upon recieving a notification, especially for recent contact, can be a response to this issue but at a high social and economic cost.\n",
    " * The \"warning after symptoms\" policy can be useful to get around the testing delays of symptomatic people if those are longer than 2/3 days. For this, though, we need to precisely identify COVID-19 symptoms and distinguish them from similar symptoms of other diseases.\n",
    " * The proportion of users of the app in the population is a central parameter. Under 20% of users, the effects of the app are insignificant, even in a best case scenario. In reality, depending on the policies adopted this threshold could be much higher.\n",
    "\n",
    "### References\n",
    "\n",
    "* [1] Contact network epidemiology: Bond percolation applied to infectious disease prediction and control, Lauren Ancel Meyers [(link)](https://www.ams.org/journals/bull/2007-44-01/S0273-0979-06-01148-7/S0273-0979-06-01148-7.pdf)\n",
    "* [2] Contact Network Epidemiology: Mathematical Methods of Modeling a Mutating Pathogen on a Two-type Network by Robert L. Seilheimer [(link)](https://repositories.lib.utexas.edu/bitstream/handle/2152/13376/Seilheimer_-_M_08.pdf)\n",
    "* [3] Network theory and SARS: predicting outbreak diversity by Lauren Ancel Meyers & al. [(link)](https://doi.org/10.1016/j.jtbi.2004.07.026)\n",
    "* [4] Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing. - PubMed - NCBI [(link)](https://science.sciencemag.org/content/sci/early/2020/04/09/science.abb6936.full.pdf)\n",
    "* [5] The Impact of Contact Tracing in Clustered Populations by Thomas House & Matt J Keeling [(link)](https://journals.plos.org/ploscompbiol/article/file?id=10.1371/journal.pcbi.1000721&type=printable)\n",
    "* [6] Temporal dynamics in viral shedding and transmissibility of COVID-19 by Xi He & al [(link)](https://www.nature.com/articles/s41591-020-0869-5)\n",
    "* [7] Clinical characteristics of 24 asymptomatic infections with COVID-19 screened among close contacts in Nanjing, China by Zhiliang Hu & al [(link)](https://link.springer.com/article/10.1007/s11427-020-1661-4)\n",
    "* [8] Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19) [(link)](https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.pdf)\n",
    "\n",
    "\n",
    "### Graph visualisation\n",
    "\n",
    "Our model is based on a graph simulation: we show a visualisation of the spread of the epidemic in the contact network.\n",
    "The contact network is organised in households (little groups of 2 to 6 nodes packed together). The members of a household have a greater probability of contact between them than with members of different households.\n",
    "\n",
    "Here is the initial state:\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/init_state.png\" alt=\"init state\" style=\"width: 500px;\"/>\n",
    "\n",
    "Color correspondence:\n",
    "\n",
    "\n",
    "| Color     | Meaning                    |\n",
    "| :-------: | :------------------------: |\n",
    "| Green     | healthy                    |\n",
    "| Red       | infected without symptoms  |\n",
    "| Pink      | infected with symptoms     |\n",
    "| Blue      | cured                      |\n",
    "| Black     | dead                       |\n",
    "\n",
    "Let's see what happens when nobody installs the app.\n",
    "\n",
    "Day 20:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/state_20_noapp.png\" alt=\"no app, day 20\" style=\"width: 500px;\"/>\n",
    "\n",
    "After 50 days, almost everyone has been infected (because almost everyone is cured):\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/state_50_no_app.png\" alt=\"no app, day 50\" style=\"width: 500px;\"/>\n",
    "\n",
    "Note that these visualisations are just meant to illustrate the model. In the \"Results\" section, we used another inital state and experimented different policies.\n",
    "\n",
    "The graph visualisations were made with Gephi and the Force Atlas 2 algorithm.\n",
    "\n",
    "### Detailed results\n",
    "\n",
    "#### Baseline case: no contact tracing\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/0%25%2C-%2C-%2C[3%2C10]%2C5j%2C0.3.png\" alt=\"no contact tracing\" style=\"width: 500px;\"/>\n",
    "\n",
    "The epidemic spreads freely from the 5% initial infected rate until almost all the population gets infected.\n",
    "The initial Rt (R0) (red curve) is a bit lower than 2 because in the initial population, 10% of the individuals are cured and because people who suffer of symptoms directly quarantine themselves.\n",
    "This was not the case at the beginning of the outbreak when R0 was evaluated around 2.\n",
    "\n",
    "We will take this graph as reference for further analysis.\n",
    "\n",
    "#### Influence of the waiting time for tests results\n",
    "\n",
    "We use these test parameters to model the current test abilities in France:\n",
    " * 5 days to wait for test results\n",
    " * a 3 to 10 days validity window\n",
    " * 30% of false positives\n",
    "\n",
    "Even with a great adoption of the app by the population (80%), no change in the way the epidemic spreads is visible because of the test delay.\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2C-%2C-%2C[3%2C10]%2C5j%2C0.3.png\" alt=\"current testing params in France, 80% adoption\" style=\"width: 500px;\"/>\n",
    "\n",
    "In fact, the critical time period to keep in mind to curb the spread of the disease is the 5 days incubation period.  \n",
    "According to the last epidemiological statistics, it seems that asymptomatic people make up to about 40% of the infected people - yet they are the cause of only 6% of infections[4].\n",
    "On the contrary, the infectiousness of presymptomatic people appears to be really high [4] [6] - so much so that most contaminations originate from them.\n",
    "Moreover, when symptoms arise the affected person goes into quarantine. This means that symptomatic transmission is contained if this measure is well respected.  \n",
    "Thus, presymptomatic contaminations are the most significant, and efforts should focus on them.\n",
    "To avoid them, it is necessary to identify and quarantine infected people at the beginning of the incubation period. Therefore, after 5-day delay it is already too late.\n",
    "\n",
    "*Note: The spikes every 5 days in test demand are due to the fact that an individual cannot ask for a new test if he is already waiting for results.*\n",
    "\n",
    "In the case of the current testing abilities in France, the characteristic time is of `2*5=10` days between the suspicions of the first infected person and the quarantining of the person notified and infected by the latter. The quarantine arrives far too late to be effective.\n",
    "\n",
    "Furthermore, the validity window of nasal tests is a key factor that limits their effectiveness. Even if we succeed in having almost perfect tests (1% of false positives) and results available with no delay, we will not be able to identify and isolate presymptomatic people that just got infected because of the 3-10 days validity window.  \n",
    "Here is a comparison with and without the validity window:\n",
    "\n",
    "* with the validity window:\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2C-%2C-%2C[3%2C10]%2C0j%2C0.01.png\" alt=\"with validity window\" style=\"width: 500px;\"/>\n",
    "* without the validity window:\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2C-%2C-%2C[-1%2C100]%2C0j%2C0.01.png\" alt=\"without validity window\" style=\"width: 500px;\"/>\n",
    "\n",
    "With the validity window, we can only identify infected people 3 days after the infection, which is already a significant part of the incubation period.\n",
    "\n",
    "* with more realistic testing parameters (2 days of delay and 15% of false negative) (still better than current tests) :\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2C-%2C-%2C[3%2C10]%2C2j%2C0.15.png\" alt=\"with better tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "Even with better screening, we can see that the curve of the test requests appears earlier, but the\n",
    "quarantine curve doesn't change much. This can be explained by the fact that a close contact has a probability of only about 2% of causing a contamination (according to a WHO report in China [8]). This leads to a high test demand from persons that only have a low chance of being infected. Moreover, in the best case, the decision to quarantine is taken 4 days after the infection, also too late.\n",
    "\n",
    "In these simulations the demand for tests is unrealistically high because we don't add any constraints on test capacity.\n",
    "We may imagine a more nuanced system to determine if we need a test: the number of notifications received, the presumed infectiousness of the infected person at the time of the contact, the comorbidity factors, etc. This could allow a more precise evaluation of the risk an individual has of getting infected.\n",
    "\n",
    "#### 'Quarantine after notification' policy\n",
    "\n",
    "* application of the \"quarantine after notification\" policy and with current tests\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2CN%2C-%2C[3%2C10]%2C5j%2C0.3.png\" alt=\"quarantine after notification with current tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "It provides a more aggressive response against the epidemic and is proven to be effective even with the current test abilities. However, a great part of its effectiveness is likely to be due to the massive quarantine it causes (80% of the population).  \n",
    "Notice the oscillations of the number of people in quarantine: when the pandemic first spreads, a great part of the population quarantines itself because of the many notifications received. After testing, a lot of them get negative results and end their quarantine, causing the epidemic to spread again.\n",
    "\n",
    "* application of \"quarantine after notification\" and with better tests (2 days of delay and 15% of false negative):\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2CN%2C-%2C[3%2C10]%2C2j%2C0.15.png\" alt=\"quarantine after notification with better tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "The effectiveness of this policy is greater with quicker tests because presymptomatic people get notified earlier.\n",
    "\n",
    "This is illustrated by the rapid collapse of Rt during the first days (red curve). \n",
    "In this case, this policy uses quick testing to get earlier and more efficient quarantines.\n",
    "With such tests parameters, the average number of quarantine days per person is also reduced as the waiting time for test results is shorter.\n",
    "\n",
    "However, it should be expected that such a policy comes with a high social cost because of its strict quarantining.\n",
    "Indeed, because of the great number of close contacts in a day (13 on average), almost all the users of the app would quarantine themselves after the first wave of contaminations.  \n",
    "This high social cost could hamper the adoption of the app.\n",
    "\n",
    "A response to this problem can be to send a notification only to latest contacts. Instead of warning all the close contacts that happended in the last 14 days when testing positive as we did in all our other simulations, we can warn only the contacts of the 7 last days. \n",
    "* application of \"quarantine after notification\" and with better tests (2 days of delay and 15% of false negative) and notificaiton of contacts in the last 7 days:\n",
    " <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%2CN%2C-%2C%5B3%2C10%5D%2C2j%2C0.15_avec_daysNotif_a_7.png\" alt=\"quarantine after notification with better tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "This way, we still warn presymptomatics in the begining of their incubation period without useless qurantine of older contacts that, even if they led to an infection, are now symptomatic or asymptomatic.\n",
    "This lowers the timer spent by healthy person in quarantine by 2 days without speeding the spread of the virus too much.\n",
    "\n",
    "In short, if test results can be available quickly, both the effectiveness of this policy and the social cost are improved.\n",
    "\n",
    "#### 'Warning after symptoms' policy\n",
    "\n",
    "* application of \"quarantine after notification\" and \"warning after symptoms\" with current tests (5 days of delay and 30% of false negatives):\n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2CN%2CW%2C[3%2C10]%2C5j%2C0.3.png\" alt=\"warning after symptoms with current tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "The proportion of notifications sent because of symptoms unrelated to COVID-19 is difficult to estimate. However if we manage to keep it fairly low ( < 0.5% chance of having COVID-looking symptoms per person per day) this policy can be an effective workaround for test delays and test reliability if applied with the \"quarantine after notification\" policy. The spread of the disease is contained at the same level as with better tests but at the cost of 50% more days of quarantine.\n",
    "\n",
    "* application of \"quarantine after notification\" and \"warning after symptoms\" with better tests (1 days of delay and 1% of false negatives): \n",
    "  <img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/80%25%2CN%2CW%2C[3%2C10]%2C1j%2C0.01.png\" alt=\"warning after symptoms with better tests\" style=\"width: 500px;\"/>\n",
    "\n",
    "In both cases it leads to the lockdown of a high proportion of the population.\n",
    "\n",
    "#### Influence of the proportion of users\n",
    "\n",
    "In the case of better tests (2 days of delay and 15% of false negatives) and the use of the \"quarantine after notification\" policy, we plot the influence of the proportion of users on several key values.\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/day_in_quarantine.png\" alt=\"days in quarantine\" style=\"width: 500px;\"/>\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/max_assymp.png\" alt=\"maximum number of asymptomatic cases\" style=\"width: 500px;\"/>\n",
    "<img src=\"https://raw.githubusercontent.com/GuilhemN/covid-19-tracing-app-modelling/master/images/proportion_of_healthy.png\" alt=\"proportion of healthy people after 60 days\" style=\"width: 500px;\"/>\n",
    "\n",
    "We see that the use of the app will not reduce the maximal pressure on healthcare system (max. symptomatic) by much but will significantly curb the total spread of the disease.\n",
    "By doubling the average number of days in quarantine per person, it triples the final proportion of healthy people.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "####################\n",
    "# GRAPH GENERATION #\n",
    "####################\n",
    "\n",
    "nbIndividuals = 1000 # number of people in the graph | nombre d'individus dans le graphe\n",
    "initHealthy = 0.85 # proportion of healthy people at start | la proportion de personnes saines à l'intant initial\n",
    "initCured = 0.1 # proportion of cured people at start | proportion de personnes guéries à l'instant initial\n",
    "# The other people are 60% presymptomatic and 40% asymptomatic at start | Les autres personnes sont 40% d'asymptomatiques et 60% de présymptomatiques au départ\n",
    "\n",
    "# graph generation for exponential degrees distribution\n",
    "#------------------------------------------------------\n",
    "deg_avg = 100 # average number of connexions per person | le nombre moyen de connexions par personne\n",
    "av_household_size = 6 # average size of household | la taille moyenne d'un foyer\n",
    "household_proba = 1 # probability of meeting a person of the same household | la probabilité de contact par jour entre membres d'un même foyer\n",
    "extern_contact_proba = 0.3 # probabilty of meeting a person of a different household | la probabilité de contact par jour entre personne de foyers différents\n",
    "\n",
    "# average contacts per day = 0.3*(100-6) + 1*6 = 34.2\n",
    "\n",
    "# graph generation with organization in households\n",
    "#-------------------------------------------------\n",
    "household_size = (3, 5) # min and max size of an household (uniform distribution) | extremums de la taille d'un foyer\n",
    "household_link = 1 # probability of contact between members of a household | proba de contact entre membres d'un foyer\n",
    "\n",
    "number_of_households = 300 # 2500 is good but a bit slow | number of households in the community | nombre de foyers dans une communauté\n",
    "community_link = 0.3 # probability of contact across households | proba de contact entre foyers\n",
    "av_deg_by_household = 400 # number of link from a household | nombre moyen de liens depuis un foyer\n",
    "\n",
    "# average external degree of an individual : 400/4 (4 is the average size of an household)\n",
    "# average contacts per day = (400/4)*0.3 + 4 = 34\n",
    "# This paramter was estimated, this is a limit of the model\n",
    "\n",
    "##############\n",
    "# APP PARAMS #\n",
    "##############\n",
    "\n",
    "daysNotif = 14 # number of days the app checks back for contact notification | nombre de jours vérifiés par l'appli pour notifier un contact\n",
    "utilApp = 0.8 # percentage of people having the app | la proportion d'utilisateurs de l'application dans la population générale\n",
    "\n",
    "pDetection = 0.9 # prob. that the app detects a contact | proba que l'appli détecte un contact\n",
    "pReport = 0.9 # prob. that a user reports his symptoms | proba qu'un utilisateur alerte de ses symptômes\n",
    "pReadNotif = 0.8 # probablity of taking a notification into account (ask for a test, quarantine) | proba de prendre en compte une notification (demande de test, quarantaine)\n",
    "\n",
    "pSymptomsNotCovid = 0.005 # every day, everyone sends a notification with prob. pSymptomsNotCovid | chaque jour, tout le monde envoie une notif avec proba PSymptomsNotCovid\n",
    "\n",
    "############\n",
    "# POLICIES #\n",
    "############\n",
    "\n",
    "# people warn the app immediately after having symptoms | on prévient l'application directement après avoir développé les symptômes\n",
    "warningAfterSymptoms = False\n",
    "\n",
    "# upon notification, an individual asks for a test (with some prob.)\n",
    "# if true, user waits for test results in quarantine, else he goes in quarantine only upon reception of positive test results\n",
    "# |\n",
    "# à la reception d'une notif, l'utilisateur demande un test (avec une certaine proba)\n",
    "# si vrai, il attend les résultats en quarantaine, sinon il ne se met en quarantaine qu'aux résultats d'un test positif\n",
    "quarantineAfterNotification = True\n",
    "\n",
    "###############\n",
    "# TEST PARAMS #\n",
    "###############\n",
    "\n",
    "testWindow = (3, 10) # tests are only effective in a given window (time since infection) | les tests ne sont efficaces que dans une fenêtre de temps après infection\n",
    "\n",
    "daysUntilResult = 2 # Time to wait before test results | attente pour l'obtention des résultats des tests\n",
    "pFalseNegative = 0.15 # prob. of false negative | proba d'avoir un faux négatif\n",
    "daysBetweenTests = 0 # the time to wait before being allowed to ask for un new test | temps avant de pouvoir demander un nouveau test\n",
    "\n",
    "\n",
    "##############\n",
    "# QUARANTINE #\n",
    "##############\n",
    "\n",
    "pQSymptoms = 0.9 # probability of going into quarantine when one has symptoms | proba de confinement lors de détection des symptômes\n",
    "quarantineFactor = 100 # reduction factor applied to the probabilities when one is in quarantine | réduction des probas de rencontre lors du confinement\n",
    "daysQuarantine = 14 # duration of the quarantine | durée de la quarantaine\n",
    "\n",
    "#################\n",
    "# PROBABILITIES #\n",
    "#################\n",
    "\n",
    "pCloseContact = 0.375 # prob. that a contact is a close contact (those detected by the app) | proba qu'un contact soit rapproché (ceux détectés par l'appli)\n",
    "pContaminationCloseContact = 0.02 # prob. of contamination after close contact with an infected person | proba de contamination après contact rapproché avec qqn d'infecté\n",
    "#according to [8] around 1 to 5% of close contact lead to virus transmission\n",
    "pContaminationCloseContactAsymp = 0.006\n",
    "# infectiousness of asymptomatic people appears to be very low according to [4] and [6]\n",
    "\n",
    "pContaminationFar = 0.001 # prob. of contamination upon non close contact (environmental or short contact) | proba de contamination par contact environnemental ou bref\n",
    "pContaminationFarAsymp = 0.0003\n",
    "\n",
    "# we took R0=2 estimate from [4] and : 34 contacts/day, an average time of infectiousness of 10 days (pre symptomatic + begining of symptoms period)\n",
    "\n",
    "#average number of infected by symptomatic : (0.375*0.02+0.625*0.001)*34*10 = 2.76\n",
    "#average number of infected by asymptomatic : (0.375*0.006+0.625*0.0003)*34*10 = 0.83\n",
    "\n",
    "# this gives 0.6*2.76 + 0.4*0.83 = 1.99 persons infected in average by an infected\n",
    "# this is plausible given the estimate of R0 and the fact that asymptomatic contamination appears to be minor\n",
    "# [4] and [6]\n",
    "\n",
    "# and (0.6*0.625*0.001 + 0.4*0.625*0.0003)*34*10 / R0 = 0.0765 -> the proportion of contaminations which are not due to close contact (environmental / short contact) (contaminations by asymptomatic people are neglected) estimated according to environmental contamination estimate in [4]\n",
    "# thus most infections (92%) are susceptible to be noticed by the app\n",
    "\n",
    "# -> the proportion of contaminations by asympt. people is : 0.4*0.83/(0.6*2.76 + 0.4*0.0.83) = 0.17 plausible according to the presumed low infectiosity shown in [4], but this is a conservative estimate (not the 0.06 given by this paper) given the high uncertainty around the results\n",
    "\n",
    "\n",
    "pAsympt = 0.4 # probability of being asymptomatic when infected | proba qu'une personne infectée soit asymptomatique\n",
    "# according to [4] and Diamond Princess estimates\n",
    "\n",
    "# parameters for the lognormal law of the incubation period | paramètres pour la loi lognormale de la période d'incubation\n",
    "incubMeanlog = 1.644 # -> ~5.5 days\n",
    "incubSdlog = 0.363 # -> ~2.1 days\n",
    "# according to [4]\n",
    "\n",
    "# The next probabilities are given for 1 step of the process, thus overall time spent in these states follows a geometric law for which expected values have been calculated\n",
    "\n",
    "pAtoG = 0.1 # probability of going from asymptomatic state to cured | proba de passer de asymptomatique à guéri\n",
    "# according to [7]\n",
    "\n",
    "pIStoC = 0.07 # probability of going from symptomatic state to cured | proba de passer de avec symptômes à gueri\n",
    "pIStoD = 0.003 # probability of dying when symptomatic | proba de décès d'une personne présentant des symptômes\n",
    "\n",
    "# average time with symptoms : 1/(0.07+0.003) = 13.7 days : plausible according to [4]\n",
    "# death rate when symptoms : 0.003/0.07 = 4.3% : plausible in France according to estimate of 1.6M cases with symptoms and 6 000 deaths the 3 April\n",
    "# https://www.mgfrance.org/publication/communiquepresse/2525-enquete-mg-france-plus-d-un-million-et-demi-de-personnes-prises-en-charge-par-leur-medecin-generaliste-pour-le-covid-19-entre-le-17-mars-et-le-3-avril\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Libs and defs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Librairies\n",
    "import random\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "# -> sliders\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "HEALTHY = 0\n",
    "ASYMP = 1\n",
    "PRESYMP = 2\n",
    "SYMP = 3\n",
    "CURED = 4\n",
    "DEAD = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Graph:\n",
    "    \"\"\" Object holding the representation of the graph and some metrics \"\"\"\n",
    "\n",
    "    def __init__(self):\n",
    "        self.individuals = []\n",
    "        self.adj = []\n",
    "\n",
    "        self.encounters = [[[] for day in range(daysNotif)] for individual in range(nbIndividuals)]\n",
    "\n",
    "        self.nbHealthy = 0 # number of healthy people\n",
    "        self.nbAS = 0 # number of asymptomatic people\n",
    "        self.nbPS = 0 # number of premptomatic people\n",
    "        self.nbS = 0 # number of symptomatic people\n",
    "        self.nbCured = 0 # number of cured persons\n",
    "        self.nbDead = 0 # number of deceased people\n",
    "        self.nbQuarantineI = 0 # number of infected people in quarantine\n",
    "        self.nbQuarantineNonI = 0 # number of non infected people in quarantine\n",
    "\n",
    "        self.nbTest = 0 # number of tests made\n",
    "\n",
    "        # cumulative counters :\n",
    "        self.nbQuarantineTotal = 0 # number of people in quarantine\n",
    "        self.nbInfectedByASPS = 0 # number of people infected by asymp. + presymp. people\n",
    "\n",
    "        #to compute Rt\n",
    "        self.stepNb = 0\n",
    "        self.contaminations = [] # number of people contaminated at a given time\n",
    "        self.numInfectedByNewInfected = [] # total number of people who will get infected by people contaminated at a given time\n",
    "\n",
    "class Individual:\n",
    "    \"\"\" Object holding the representation of an individual \"\"\"\n",
    "\n",
    "    def __init__(self, state, daysQuarantine, app, sentNotification, daysIncubation, timeSinceInfection, timeLeftForTestResult):\n",
    "        self.state = state\n",
    "        self.daysQuarantine = daysQuarantine\n",
    "        self.app = app\n",
    "        self.sentNotification = sentNotification\n",
    "        self.daysIncubation = daysIncubation\n",
    "        self.timeSinceInfection = timeSinceInfection\n",
    "        self.timeSinceLastTest = np.inf # we don't want to test people too often\n",
    "        self.timeLeftForTestResult = timeLeftForTestResult\n",
    "        self.nbInfected = 0\n",
    "\n",
    "    def in_state(self, state):\n",
    "        return self.state == state\n",
    "\n",
    "    def is_infected(self):\n",
    "        return self.state in [PRESYMP, ASYMP, SYMP]\n",
    "\n",
    "    def has_no_covid(self):\n",
    "        return self.state in [HEALTHY, CURED]\n",
    "\n",
    "    def in_quarantine(self):\n",
    "        return self.daysQuarantine > 0\n",
    "\n",
    "    def go_quarantine(self):\n",
    "        if self.daysQuarantine <= 0:\n",
    "            self.daysQuarantine = daysQuarantine # goes into quarantine if isn't already"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def create_individuals(graph):\n",
    "    graph.contaminations.append(0)\n",
    "    for i in range(nbIndividuals):\n",
    "        app = False\n",
    "        if random.uniform(0,1) < utilApp:\n",
    "            app = True\n",
    "        s = PRESYMP\n",
    "        time_since_infection = -1\n",
    "        incub = 0\n",
    "        r = random.random()\n",
    "        if r < initHealthy:\n",
    "            s = HEALTHY\n",
    "            graph.nbHealthy += 1\n",
    "        elif r < initHealthy + initCured:\n",
    "            s = CURED\n",
    "            graph.nbCured += 1\n",
    "        else:\n",
    "            graph.contaminations[0] += 1 # we start as if a proportion of the population just got infected\n",
    "            time_since_infection = 0\n",
    "            if random.random() < pAsympt:\n",
    "                s = ASYMP\n",
    "                graph.nbAS += 1\n",
    "            else:\n",
    "                s = PRESYMP\n",
    "                incub = round(np.random.lognormal(incubMeanlog, incubSdlog))\n",
    "                graph.nbPS += 1\n",
    "\n",
    "        # state, quarantine, app, notif, incubation, timeSinceInfection, timeLeftForTestResult\n",
    "        graph.individuals.append(Individual(s,  0, app, False, incub, time_since_infection, -1))\n",
    "\n",
    "\n",
    "def init_graph_exp(graph):\n",
    "    \"\"\" Graph initialisation based on exponential ditribution of degrees \"\"\"\n",
    "\n",
    "    create_individuals(graph)\n",
    "\n",
    "    # affecting degrees to vertices\n",
    "    degrees = np.around(np.random.exponential(deg_avg, nbIndividuals))\n",
    "\n",
    "    # to get an even number of total degrees\n",
    "    S = sum(degrees)\n",
    "    if S%2 == 1:\n",
    "        degrees[0] += 1\n",
    "        S += 1\n",
    "\n",
    "    graph.adj = [[] for i in range(nbIndividuals)]\n",
    "    while S > 0:\n",
    "        # creating an edge\n",
    "        [p1, p2] = np.random.choice(len(degrees), 2, replace=False, p=degrees/S)\n",
    "        if degrees[p1] <= av_household_size or degrees[p2] <= av_household_size:\n",
    "        \t# the last edges created are edges within households\n",
    "            graph.adj[p1].append({\"node\" : p2, \"proba\" : household_proba})\n",
    "            graph.adj[p2].append({\"node\" : p1, \"proba\" : household_proba})\n",
    "        else:\n",
    "            graph.adj[p1].append({\"node\" : p2, \"proba\" : extern_contact_proba})\n",
    "            graph.adj[p2].append({\"node\" : p1, \"proba\" : extern_contact_proba})\n",
    "        degrees[p1] -= 1\n",
    "        degrees[p2] -= 1\n",
    "        S -= 2\n",
    "\n",
    "\n",
    "def init_graph_household(graph):\n",
    "    \"\"\" Graph generation based on households organisation \"\"\"\n",
    "\n",
    "    global nbIndividuals\n",
    "\n",
    "    # creation of the households\n",
    "    graph.adj = []\n",
    "\n",
    "    for i in range(number_of_households):\n",
    "        size = random.randint(household_size[0], household_size[1])\n",
    "        nb = len(graph.adj)\n",
    "        for i in range(nb, nb+size):\n",
    "            household = []\n",
    "            for j in range(nb, nb+size):\n",
    "                if (i != j):\n",
    "                    household.append({\"node\": j, \"proba\": household_link})\n",
    "            graph.adj.append(household)\n",
    "\n",
    "    # linkage of the households\n",
    "    for i in range(av_deg_by_household*number_of_households):\n",
    "        [p1, p2] = np.random.choice(len(graph.adj), 2, replace=False)\n",
    "        graph.adj[p1].append({\"node\": p2, \"proba\": community_link})\n",
    "        graph.adj[p2].append({\"node\": p1, \"proba\": community_link})\n",
    "\n",
    "    nbIndividuals = len(graph.adj)\n",
    "\n",
    "    create_individuals(graph)\n",
    "\n",
    "    graph.encounters = [[[] for day in range(daysNotif)] for individual in range(nbIndividuals)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Updating the graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contamination(graph, i, j, closeContact):\n",
    "    \"\"\" Individuals i and j have come into contact, leading to a possible contamination | Les individus i et j sont entrés en contact, une contamination est possible \"\"\"\n",
    "\n",
    "    if graph.individuals[i].state == graph.individuals[j].state:\n",
    "        return\n",
    "\n",
    "    if graph.individuals[i].in_state(HEALTHY):\n",
    "        contamination(graph, j, i, closeContact)\n",
    "        return\n",
    "\n",
    "    # i is the infected individual\n",
    "    if graph.individuals[i].is_infected():\n",
    "        if graph.individuals[j].in_state(HEALTHY):\n",
    "\n",
    "            if closeContact:\n",
    "                pContamination = pContaminationCloseContact\n",
    "                pContaminationAsymp = pContaminationCloseContactAsymp\n",
    "            else:\n",
    "                pContamination = pContaminationFar\n",
    "                pContaminationAsymp = pContaminationFarAsymp\n",
    "\n",
    "            if (random.random() < pContamination and (not graph.individuals[i].in_state(ASYMP))) or \\\n",
    "                (random.random() < pContaminationAsymp and graph.individuals[i].in_state(ASYMP)):\n",
    "                # j becomes infected\n",
    "\n",
    "                # for Rt computation\n",
    "                graph.contaminations[graph.stepNb] += 1\n",
    "                graph.numInfectedByNewInfected[graph.stepNb - graph.individuals[i].timeSinceInfection] += 1 # parent infection took place timeSinceInfection ago\n",
    "\n",
    "                if graph.individuals[i].in_state(ASYMP) or graph.individuals[i].in_state(PRESYMP):\n",
    "                    graph.nbInfectedByASPS += 1\n",
    "                graph.individuals[j].timeSinceInfection = 0\n",
    "                graph.individuals[i].nbInfected += 1 # i has infected one more person\n",
    "                graph.nbHealthy -= 1\n",
    "                if random.random() < pAsympt:\n",
    "                    graph.individuals[j].state = ASYMP\n",
    "                    graph.nbAS += 1\n",
    "                else:\n",
    "                    graph.individuals[j].state = PRESYMP\n",
    "                    graph.individuals[j].daysIncubation = round(np.random.lognormal(incubMeanlog, incubSdlog))\n",
    "                    graph.nbPS += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def test_individual(individual, graph):\n",
    "    # if there is a test incoming, the person is not tested again\n",
    "    if individual.timeLeftForTestResult >= 0 or individual.in_state(DEAD):\n",
    "        return\n",
    "\n",
    "    # the person was tested not long ago\n",
    "    if individual.timeSinceLastTest < daysBetweenTests:\n",
    "        return\n",
    "\n",
    "    # the person is tested\n",
    "    individual.timeSinceLastTest = 0\n",
    "    graph.nbTest += 1\n",
    "    individual.timeLeftForTestResult = daysUntilResult\n",
    "\n",
    "    if individual.has_no_covid():\n",
    "        individual.latestTestResult = False # we assume that there are no false positives\n",
    "        return\n",
    "\n",
    "    if individual.timeSinceInfection < testWindow[0] or individual.timeSinceInfection > testWindow[1]:\n",
    "        individual.latestTestResult = False # not in the detection window, the test fails\n",
    "        return\n",
    "\n",
    "    # otherwise the person is ill\n",
    "    # the test result depends whether we have a false negative or not\n",
    "    individual.latestTestResult = not (random.random() < pFalseNegative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def send_notification(graph, i):\n",
    "    \"\"\" Send notification to people who have been in touch with i | Envoi d'une notif aux personnes ayant été en contact avec i \"\"\"\n",
    "\n",
    "    if graph.individuals[i].sentNotification:\n",
    "        return # notifications already sent\n",
    "\n",
    "    graph.individuals[i].sentNotification = True\n",
    "    for daysEncounter in graph.encounters[i]:\n",
    "        # note: graph.encounter[i] is empty if i does not have the app so there is no need to have an additional condition\n",
    "        for contact in daysEncounter:\n",
    "            if random.random() < pReadNotif: # if the person takes the notification into account\n",
    "                # the person is always tested \n",
    "                test_individual(graph.individuals[contact], graph) # asks for a test\n",
    "                if quarantineAfterNotification: # in this case, the person waits for test results in quarantine\n",
    "                    graph.individuals[contact].go_quarantine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def make_encounters(graph, i):\n",
    "    \"\"\" Assess all encounters made by i in one day | Détermine toutes les rencontres faites par i en un jour \"\"\"\n",
    "\n",
    "    for edge in graph.adj[i]:\n",
    "        j = edge['node']\n",
    "        if j < i:\n",
    "            continue # only check one way of the edge | on ne regarde qu'un sens de chaque arête\n",
    "\n",
    "        # if i and/or j are in quarantine, reduce the probability that they meet | si i et/ou j sont confinés, réduction de leur proba de rencontre\n",
    "        factor = 1\n",
    "        if graph.individuals[i].in_quarantine():\n",
    "            factor *= quarantineFactor\n",
    "        if graph.individuals[j].in_quarantine():\n",
    "            factor *= quarantineFactor\n",
    "\n",
    "        if random.random() < edge['proba'] / factor:\n",
    "            if random.random() < pCloseContact: # if this is a close contact\n",
    "                # if i and j have the app, we save their encounter | si i et j ont l'appli, on note la rencontre\n",
    "                if graph.individuals[i].app and graph.individuals[j].app and random.random() < pDetection: # contact detections are symmetric in our model\n",
    "                    graph.encounters[i][-1].append(j)\n",
    "                    graph.encounters[j][-1].append(i)\n",
    "                contamination(graph, i, j, True)\n",
    "            else:\n",
    "                contamination(graph, i, j, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def step(graph):\n",
    "    \"\"\" Step from a day to the next day | Passage au jour suivant du graphe \"\"\"\n",
    "\n",
    "    graph.nbTest = 0\n",
    "    for encounter in graph.encounters:\n",
    "        encounter.append([]) # will contain every encounter of the day | contiendra les nouvelles rencontres du jour\n",
    "\n",
    "    graph.contaminations.append(0)\n",
    "    graph.numInfectedByNewInfected.append(0)\n",
    "\n",
    "    ## go through each possible encounter | on constate toutes les rencontres entre individus\n",
    "    for i in range(nbIndividuals):\n",
    "        make_encounters(graph, i)\n",
    "\n",
    "    ## update the states | on met à jour les états des individus\n",
    "    for i, individual in enumerate(graph.individuals):\n",
    "        if individual.in_state(ASYMP):\n",
    "            if random.random() < pAtoG:\n",
    "                graph.nbAS -= 1\n",
    "                graph.nbCured += 1\n",
    "                individual.state = CURED\n",
    "\n",
    "        elif individual.in_state(PRESYMP):\n",
    "            if individual.daysIncubation == 0: # the person develops symptoms\n",
    "                graph.nbPS -= 1\n",
    "                graph.nbS += 1\n",
    "                individual.state = SYMP\n",
    "\n",
    "                # send the notifications (encounters[i] is empty if i doesn't have the app) | envoi des notifs (encounters[i] vide si i n'a pas l'appli)\n",
    "                if random.random() < pReport and warningAfterSymptoms:\n",
    "                    send_notification(graph, i)\n",
    "                if random.random() < pQSymptoms: # go into quarantine if symptoms appear | mise en confinement à la détection des symptômes\n",
    "                    individual.go_quarantine()\n",
    "\n",
    "                test_individual(individual, graph) # all individuals developing symptoms are tested \n",
    "\n",
    "        elif individual.in_state(SYMP):\n",
    "            action = random.random()\n",
    "            if action < pIStoC:\n",
    "                graph.nbS -= 1\n",
    "                graph.nbCured += 1\n",
    "                individual.state = CURED\n",
    "            elif action > 1 - pIStoD:\n",
    "                graph.nbS -= 1\n",
    "                graph.nbDead += 1\n",
    "                individual.state = DEAD\n",
    "\n",
    "        # if warningAfterSymptoms is True, each individual has a probability of sending a false notification due to symptoms that are misinterpreted as from COVID-19\n",
    "        # | si warningAfterSymptoms est vrai, chaque individu a une probabilité d'envoyer une notification en raison de symptômes faussement perçus comme relevant du COVID-19\n",
    "        if warningAfterSymptoms and random.random() < pSymptomsNotCovid:\n",
    "            send_notification(graph, i)\n",
    "\n",
    "        # reception of test results | réception des résultats de test\n",
    "        if individual.timeLeftForTestResult == 0:\n",
    "            if individual.in_quarantine() and individual.latestTestResult == False: # is in quarantine and gets a negative test\n",
    "                individual.daysQuarantine = 0 # end of quarantine\n",
    "\n",
    "            if individual.latestTestResult == True:\n",
    "                individual.go_quarantine()\n",
    "                individual.timeLeftForTestResult = np.inf # people tested positive are not tested again\n",
    "\n",
    "                if random.random() < pReport: # not everyone reports a positive test to the app\n",
    "                    send_notification(graph, i)\n",
    "\n",
    "                individual.app = False # unsubscribe from the app in order to not consider new notifications\n",
    "        individual.timeLeftForTestResult -= 1\n",
    "\n",
    "    ## results of the day | bilan du jour\n",
    "    graph.nbQuarantineNonI = 0\n",
    "    graph.nbQuarantineI = 0\n",
    "\n",
    "    for individual in graph.individuals:\n",
    "        if individual.in_state(DEAD):\n",
    "            continue\n",
    "\n",
    "        individual.daysQuarantine -= 1\n",
    "        individual.daysIncubation -= 1\n",
    "        individual.timeSinceLastTest += 1\n",
    "\n",
    "        # if there are still symptoms we don't end the quarantine\n",
    "        if (not individual.in_quarantine()) and individual.in_state(SYMP):\n",
    "            individual.daysQuarantine = 1\n",
    "\n",
    "        if individual.in_quarantine():\n",
    "            graph.nbQuarantineTotal += 1/nbIndividuals\n",
    "            if not individual.is_infected():\n",
    "                graph.nbQuarantineNonI += 1\n",
    "            else:\n",
    "                graph.nbQuarantineI += 1\n",
    "\n",
    "        if individual.timeSinceInfection >= 0:\n",
    "            individual.timeSinceInfection += 1\n",
    "\n",
    "    ## deleting oldest recorded day | suppression du plus vieux jour de l'historique\n",
    "    for encounter in graph.encounters:\n",
    "        encounter.pop(0)\n",
    "    graph.stepNb += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Display\n",
    "Interactive model below (it takes about 10-15 sec to appear and to run a simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<IPython.core.display.Javascript object>",
      "application/javascript": "/* Put everything inside the global mpl namespace */\nwindow.mpl = {};\n\n\nmpl.get_websocket_type = function() {\n    if (typeof(WebSocket) !== 'undefined') {\n        return WebSocket;\n    } else if (typeof(MozWebSocket) !== 'undefined') {\n        return MozWebSocket;\n    } else {\n        alert('Your browser does not have WebSocket support. ' +\n              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n              'Firefox 4 and 5 are also supported but you ' +\n              'have to enable WebSockets in about:config.');\n    };\n}\n\nmpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n    this.id = figure_id;\n\n    this.ws = websocket;\n\n    this.supports_binary = (this.ws.binaryType != undefined);\n\n    if (!this.supports_binary) {\n        var warnings = document.getElementById(\"mpl-warnings\");\n        if (warnings) {\n            warnings.style.display = 'block';\n            warnings.textContent = (\n                \"This browser does not support binary websocket messages. \" +\n                    \"Performance may be slow.\");\n        }\n    }\n\n    this.imageObj = new Image();\n\n    this.context = undefined;\n    this.message = undefined;\n    this.canvas = undefined;\n    this.rubberband_canvas = undefined;\n    this.rubberband_context = undefined;\n    this.format_dropdown = undefined;\n\n    this.image_mode = 'full';\n\n    this.root = $('<div/>');\n    this._root_extra_style(this.root)\n    this.root.attr('style', 'display: inline-block');\n\n    $(parent_element).append(this.root);\n\n    this._init_header(this);\n    this._init_canvas(this);\n    this._init_toolbar(this);\n\n    var fig = this;\n\n    this.waiting = false;\n\n    this.ws.onopen =  function () {\n            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n            fig.send_message(\"send_image_mode\", {});\n            if (mpl.ratio != 1) {\n                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n            }\n            fig.send_message(\"refresh\", {});\n        }\n\n    this.imageObj.onload = function() {\n            if (fig.image_mode == 'full') {\n                // Full images could contain transparency (where diff images\n                // almost always do), so we need to clear the canvas so that\n                // there is no ghosting.\n                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n            }\n            fig.context.drawImage(fig.imageObj, 0, 0);\n        };\n\n    this.imageObj.onunload = function() {\n        fig.ws.close();\n    }\n\n    this.ws.onmessage = this._make_on_message_function(this);\n\n    this.ondownload = ondownload;\n}\n\nmpl.figure.prototype._init_header = function() {\n    var titlebar = $(\n        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n        'ui-helper-clearfix\"/>');\n    var titletext = $(\n        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n        'text-align: center; padding: 3px;\"/>');\n    titlebar.append(titletext)\n    this.root.append(titlebar);\n    this.header = titletext[0];\n}\n\n\n\nmpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n\n}\n\n\nmpl.figure.prototype._root_extra_style = function(canvas_div) {\n\n}\n\nmpl.figure.prototype._init_canvas = function() {\n    var fig = this;\n\n    var canvas_div = $('<div/>');\n\n    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n\n    function canvas_keyboard_event(event) {\n        return fig.key_event(event, event['data']);\n    }\n\n    canvas_div.keydown('key_press', canvas_keyboard_event);\n    canvas_div.keyup('key_release', canvas_keyboard_event);\n    this.canvas_div = canvas_div\n    this._canvas_extra_style(canvas_div)\n    this.root.append(canvas_div);\n\n    var canvas = $('<canvas/>');\n    canvas.addClass('mpl-canvas');\n    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n\n    this.canvas = canvas[0];\n    this.context = canvas[0].getContext(\"2d\");\n\n    var backingStore = this.context.backingStorePixelRatio ||\n\tthis.context.webkitBackingStorePixelRatio ||\n\tthis.context.mozBackingStorePixelRatio ||\n\tthis.context.msBackingStorePixelRatio ||\n\tthis.context.oBackingStorePixelRatio ||\n\tthis.context.backingStorePixelRatio || 1;\n\n    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n\n    var rubberband = $('<canvas/>');\n    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n\n    var pass_mouse_events = true;\n\n    canvas_div.resizable({\n        start: function(event, ui) {\n            pass_mouse_events = false;\n        },\n        resize: function(event, ui) {\n            fig.request_resize(ui.size.width, ui.size.height);\n        },\n        stop: function(event, ui) {\n            pass_mouse_events = true;\n            fig.request_resize(ui.size.width, ui.size.height);\n        },\n    });\n\n    function mouse_event_fn(event) {\n        if (pass_mouse_events)\n            return fig.mouse_event(event, event['data']);\n    }\n\n    rubberband.mousedown('button_press', mouse_event_fn);\n    rubberband.mouseup('button_release', mouse_event_fn);\n    // Throttle sequential mouse events to 1 every 20ms.\n    rubberband.mousemove('motion_notify', mouse_event_fn);\n\n    rubberband.mouseenter('figure_enter', mouse_event_fn);\n    rubberband.mouseleave('figure_leave', mouse_event_fn);\n\n    canvas_div.on(\"wheel\", function (event) {\n        event = event.originalEvent;\n        event['data'] = 'scroll'\n        if (event.deltaY < 0) {\n            event.step = 1;\n        } else {\n            event.step = -1;\n        }\n        mouse_event_fn(event);\n    });\n\n    canvas_div.append(canvas);\n    canvas_div.append(rubberband);\n\n    this.rubberband = rubberband;\n    this.rubberband_canvas = rubberband[0];\n    this.rubberband_context = rubberband[0].getContext(\"2d\");\n    this.rubberband_context.strokeStyle = \"#000000\";\n\n    this._resize_canvas = function(width, height) {\n        // Keep the size of the canvas, canvas container, and rubber band\n        // canvas in synch.\n        canvas_div.css('width', width)\n        canvas_div.css('height', height)\n\n        canvas.attr('width', width * mpl.ratio);\n        canvas.attr('height', height * mpl.ratio);\n        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n\n        rubberband.attr('width', width);\n        rubberband.attr('height', height);\n    }\n\n    // Set the figure to an initial 600x600px, this will subsequently be updated\n    // upon first draw.\n    this._resize_canvas(600, 600);\n\n    // Disable right mouse context menu.\n    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n        return false;\n    });\n\n    function set_focus () {\n        canvas.focus();\n        canvas_div.focus();\n    }\n\n    window.setTimeout(set_focus, 100);\n}\n\nmpl.figure.prototype._init_toolbar = function() {\n    var fig = this;\n\n    var nav_element = $('<div/>');\n    nav_element.attr('style', 'width: 100%');\n    this.root.append(nav_element);\n\n    // Define a callback function for later on.\n    function toolbar_event(event) {\n        return fig.toolbar_button_onclick(event['data']);\n    }\n    function toolbar_mouse_event(event) {\n        return fig.toolbar_button_onmouseover(event['data']);\n    }\n\n    for(var toolbar_ind in mpl.toolbar_items) {\n        var name = mpl.toolbar_items[toolbar_ind][0];\n        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n        var image = mpl.toolbar_items[toolbar_ind][2];\n        var method_name = mpl.toolbar_items[toolbar_ind][3];\n\n        if (!name) {\n            // put a spacer in here.\n            continue;\n        }\n        var button = $('<button/>');\n        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n                        'ui-button-icon-only');\n        button.attr('role', 'button');\n        button.attr('aria-disabled', 'false');\n        button.click(method_name, toolbar_event);\n        button.mouseover(tooltip, toolbar_mouse_event);\n\n        var icon_img = $('<span/>');\n        icon_img.addClass('ui-button-icon-primary ui-icon');\n        icon_img.addClass(image);\n        icon_img.addClass('ui-corner-all');\n\n        var tooltip_span = $('<span/>');\n        tooltip_span.addClass('ui-button-text');\n        tooltip_span.html(tooltip);\n\n        button.append(icon_img);\n        button.append(tooltip_span);\n\n        nav_element.append(button);\n    }\n\n    var fmt_picker_span = $('<span/>');\n\n    var fmt_picker = $('<select/>');\n    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n    fmt_picker_span.append(fmt_picker);\n    nav_element.append(fmt_picker_span);\n    this.format_dropdown = fmt_picker[0];\n\n    for (var ind in mpl.extensions) {\n        var fmt = mpl.extensions[ind];\n        var option = $(\n            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n        fmt_picker.append(option);\n    }\n\n    // Add hover states to the ui-buttons\n    $( \".ui-button\" ).hover(\n        function() { $(this).addClass(\"ui-state-hover\");},\n        function() { $(this).removeClass(\"ui-state-hover\");}\n    );\n\n    var status_bar = $('<span class=\"mpl-message\"/>');\n    nav_element.append(status_bar);\n    this.message = status_bar[0];\n}\n\nmpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n    // which will in turn request a refresh of the image.\n    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n}\n\nmpl.figure.prototype.send_message = function(type, properties) {\n    properties['type'] = type;\n    properties['figure_id'] = this.id;\n    this.ws.send(JSON.stringify(properties));\n}\n\nmpl.figure.prototype.send_draw_message = function() {\n    if (!this.waiting) {\n        this.waiting = true;\n        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n    }\n}\n\n\nmpl.figure.prototype.handle_save = function(fig, msg) {\n    var format_dropdown = fig.format_dropdown;\n    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n    fig.ondownload(fig, format);\n}\n\n\nmpl.figure.prototype.handle_resize = function(fig, msg) {\n    var size = msg['size'];\n    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n        fig._resize_canvas(size[0], size[1]);\n        fig.send_message(\"refresh\", {});\n    };\n}\n\nmpl.figure.prototype.handle_rubberband = function(fig, msg) {\n    var x0 = msg['x0'] / mpl.ratio;\n    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n    var x1 = msg['x1'] / mpl.ratio;\n    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n    x0 = Math.floor(x0) + 0.5;\n    y0 = Math.floor(y0) + 0.5;\n    x1 = Math.floor(x1) + 0.5;\n    y1 = Math.floor(y1) + 0.5;\n    var min_x = Math.min(x0, x1);\n    var min_y = Math.min(y0, y1);\n    var width = Math.abs(x1 - x0);\n    var height = Math.abs(y1 - y0);\n\n    fig.rubberband_context.clearRect(\n        0, 0, fig.canvas.width / mpl.ratio, fig.canvas.height / mpl.ratio);\n\n    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n}\n\nmpl.figure.prototype.handle_figure_label = function(fig, msg) {\n    // Updates the figure title.\n    fig.header.textContent = msg['label'];\n}\n\nmpl.figure.prototype.handle_cursor = function(fig, msg) {\n    var cursor = msg['cursor'];\n    switch(cursor)\n    {\n    case 0:\n        cursor = 'pointer';\n        break;\n    case 1:\n        cursor = 'default';\n        break;\n    case 2:\n        cursor = 'crosshair';\n        break;\n    case 3:\n        cursor = 'move';\n        break;\n    }\n    fig.rubberband_canvas.style.cursor = cursor;\n}\n\nmpl.figure.prototype.handle_message = function(fig, msg) {\n    fig.message.textContent = msg['message'];\n}\n\nmpl.figure.prototype.handle_draw = function(fig, msg) {\n    // Request the server to send over a new figure.\n    fig.send_draw_message();\n}\n\nmpl.figure.prototype.handle_image_mode = function(fig, msg) {\n    fig.image_mode = msg['mode'];\n}\n\nmpl.figure.prototype.updated_canvas_event = function() {\n    // Called whenever the canvas gets updated.\n    this.send_message(\"ack\", {});\n}\n\n// A function to construct a web socket function for onmessage handling.\n// Called in the figure constructor.\nmpl.figure.prototype._make_on_message_function = function(fig) {\n    return function socket_on_message(evt) {\n        if (evt.data instanceof Blob) {\n            /* FIXME: We get \"Resource interpreted as Image but\n             * transferred with MIME type text/plain:\" errors on\n             * Chrome.  But how to set the MIME type?  It doesn't seem\n             * to be part of the websocket stream */\n            evt.data.type = \"image/png\";\n\n            /* Free the memory for the previous frames */\n            if (fig.imageObj.src) {\n                (window.URL || window.webkitURL).revokeObjectURL(\n                    fig.imageObj.src);\n            }\n\n            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n                evt.data);\n            fig.updated_canvas_event();\n            fig.waiting = false;\n            return;\n        }\n        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n            fig.imageObj.src = evt.data;\n            fig.updated_canvas_event();\n            fig.waiting = false;\n            return;\n        }\n\n        var msg = JSON.parse(evt.data);\n        var msg_type = msg['type'];\n\n        // Call the  \"handle_{type}\" callback, which takes\n        // the figure and JSON message as its only arguments.\n        try {\n            var callback = fig[\"handle_\" + msg_type];\n        } catch (e) {\n            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n            return;\n        }\n\n        if (callback) {\n            try {\n                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n                callback(fig, msg);\n            } catch (e) {\n                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n            }\n        }\n    };\n}\n\n// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\nmpl.findpos = function(e) {\n    //this section is from http://www.quirksmode.org/js/events_properties.html\n    var targ;\n    if (!e)\n        e = window.event;\n    if (e.target)\n        targ = e.target;\n    else if (e.srcElement)\n        targ = e.srcElement;\n    if (targ.nodeType == 3) // defeat Safari bug\n        targ = targ.parentNode;\n\n    // jQuery normalizes the pageX and pageY\n    // pageX,Y are the mouse positions relative to the document\n    // offset() returns the position of the element relative to the document\n    var x = e.pageX - $(targ).offset().left;\n    var y = e.pageY - $(targ).offset().top;\n\n    return {\"x\": x, \"y\": y};\n};\n\n/*\n * return a copy of an object with only non-object keys\n * we need this to avoid circular references\n * http://stackoverflow.com/a/24161582/3208463\n */\nfunction simpleKeys (original) {\n  return Object.keys(original).reduce(function (obj, key) {\n    if (typeof original[key] !== 'object')\n        obj[key] = original[key]\n    return obj;\n  }, {});\n}\n\nmpl.figure.prototype.mouse_event = function(event, name) {\n    var canvas_pos = mpl.findpos(event)\n\n    if (name === 'button_press')\n    {\n        this.canvas.focus();\n        this.canvas_div.focus();\n    }\n\n    var x = canvas_pos.x * mpl.ratio;\n    var y = canvas_pos.y * mpl.ratio;\n\n    this.send_message(name, {x: x, y: y, button: event.button,\n                             step: event.step,\n                             guiEvent: simpleKeys(event)});\n\n    /* This prevents the web browser from automatically changing to\n     * the text insertion cursor when the button is pressed.  We want\n     * to control all of the cursor setting manually through the\n     * 'cursor' event from matplotlib */\n    event.preventDefault();\n    return false;\n}\n\nmpl.figure.prototype._key_event_extra = function(event, name) {\n    // Handle any extra behaviour associated with a key event\n}\n\nmpl.figure.prototype.key_event = function(event, name) {\n\n    // Prevent repeat events\n    if (name == 'key_press')\n    {\n        if (event.which === this._key)\n            return;\n        else\n            this._key = event.which;\n    }\n    if (name == 'key_release')\n        this._key = null;\n\n    var value = '';\n    if (event.ctrlKey && event.which != 17)\n        value += \"ctrl+\";\n    if (event.altKey && event.which != 18)\n        value += \"alt+\";\n    if (event.shiftKey && event.which != 16)\n        value += \"shift+\";\n\n    value += 'k';\n    value += event.which.toString();\n\n    this._key_event_extra(event, name);\n\n    this.send_message(name, {key: value,\n                             guiEvent: simpleKeys(event)});\n    return false;\n}\n\nmpl.figure.prototype.toolbar_button_onclick = function(name) {\n    if (name == 'download') {\n        this.handle_save(this, null);\n    } else {\n        this.send_message(\"toolbar_button\", {name: name});\n    }\n};\n\nmpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n    this.message.textContent = tooltip;\n};\nmpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n\nmpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n\nmpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n    // Create a \"websocket\"-like object which calls the given IPython comm\n    // object with the appropriate methods. Currently this is a non binary\n    // socket, so there is still some room for performance tuning.\n    var ws = {};\n\n    ws.close = function() {\n        comm.close()\n    };\n    ws.send = function(m) {\n        //console.log('sending', m);\n        comm.send(m);\n    };\n    // Register the callback with on_msg.\n    comm.on_msg(function(msg) {\n        //console.log('receiving', msg['content']['data'], msg);\n        // Pass the mpl event to the overridden (by mpl) onmessage function.\n        ws.onmessage(msg['content']['data'])\n    });\n    return ws;\n}\n\nmpl.mpl_figure_comm = function(comm, msg) {\n    // This is the function which gets called when the mpl process\n    // starts-up an IPython Comm through the \"matplotlib\" channel.\n\n    var id = msg.content.data.id;\n    // Get hold of the div created by the display call when the Comm\n    // socket was opened in Python.\n    var element = $(\"#\" + id);\n    var ws_proxy = comm_websocket_adapter(comm)\n\n    function ondownload(figure, format) {\n        window.open(figure.imageObj.src);\n    }\n\n    var fig = new mpl.figure(id, ws_proxy,\n                           ondownload,\n                           element.get(0));\n\n    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n    // web socket which is closed, not our websocket->open comm proxy.\n    ws_proxy.onopen();\n\n    fig.parent_element = element.get(0);\n    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n    if (!fig.cell_info) {\n        console.error(\"Failed to find cell for figure\", id, fig);\n        return;\n    }\n\n    var output_index = fig.cell_info[2]\n    var cell = fig.cell_info[0];\n\n};\n\nmpl.figure.prototype.handle_close = function(fig, msg) {\n    var width = fig.canvas.width/mpl.ratio\n    fig.root.unbind('remove')\n\n    // Update the output cell to use the data from the current canvas.\n    fig.push_to_output();\n    var dataURL = fig.canvas.toDataURL();\n    // Re-enable the keyboard manager in IPython - without this line, in FF,\n    // the notebook keyboard shortcuts fail.\n    IPython.keyboard_manager.enable()\n    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n    fig.close_ws(fig, msg);\n}\n\nmpl.figure.prototype.close_ws = function(fig, msg){\n    fig.send_message('closing', msg);\n    // fig.ws.close()\n}\n\nmpl.figure.prototype.push_to_output = function(remove_interactive) {\n    // Turn the data on the canvas into data in the output cell.\n    var width = this.canvas.width/mpl.ratio\n    var dataURL = this.canvas.toDataURL();\n    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n}\n\nmpl.figure.prototype.updated_canvas_event = function() {\n    // Tell IPython that the notebook contents must change.\n    IPython.notebook.set_dirty(true);\n    this.send_message(\"ack\", {});\n    var fig = this;\n    // Wait a second, then push the new image to the DOM so\n    // that it is saved nicely (might be nice to debounce this).\n    setTimeout(function () { fig.push_to_output() }, 1000);\n}\n\nmpl.figure.prototype._init_toolbar = function() {\n    var fig = this;\n\n    var nav_element = $('<div/>');\n    nav_element.attr('style', 'width: 100%');\n    this.root.append(nav_element);\n\n    // Define a callback function for later on.\n    function toolbar_event(event) {\n        return fig.toolbar_button_onclick(event['data']);\n    }\n    function toolbar_mouse_event(event) {\n        return fig.toolbar_button_onmouseover(event['data']);\n    }\n\n    for(var toolbar_ind in mpl.toolbar_items){\n        var name = mpl.toolbar_items[toolbar_ind][0];\n        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n        var image = mpl.toolbar_items[toolbar_ind][2];\n        var method_name = mpl.toolbar_items[toolbar_ind][3];\n\n        if (!name) { continue; };\n\n        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n        button.click(method_name, toolbar_event);\n        button.mouseover(tooltip, toolbar_mouse_event);\n        nav_element.append(button);\n    }\n\n    // Add the status bar.\n    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n    nav_element.append(status_bar);\n    this.message = status_bar[0];\n\n    // Add the close button to the window.\n    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n    button.click(function (evt) { fig.handle_close(fig, {}); } );\n    button.mouseover('Stop Interaction', toolbar_mouse_event);\n    buttongrp.append(button);\n    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n    titlebar.prepend(buttongrp);\n}\n\nmpl.figure.prototype._root_extra_style = function(el){\n    var fig = this\n    el.on(\"remove\", function(){\n\tfig.close_ws(fig, {});\n    });\n}\n\nmpl.figure.prototype._canvas_extra_style = function(el){\n    // this is important to make the div 'focusable\n    el.attr('tabindex', 0)\n    // reach out to IPython and tell the keyboard manager to turn it's self\n    // off when our div gets focus\n\n    // location in version 3\n    if (IPython.notebook.keyboard_manager) {\n        IPython.notebook.keyboard_manager.register_events(el);\n    }\n    else {\n        // location in version 2\n        IPython.keyboard_manager.register_events(el);\n    }\n\n}\n\nmpl.figure.prototype._key_event_extra = function(event, name) {\n    var manager = IPython.notebook.keyboard_manager;\n    if (!manager)\n        manager = IPython.keyboard_manager;\n\n    // Check for shift+enter\n    if (event.shiftKey && event.which == 13) {\n        this.canvas_div.blur();\n        // select the cell after this one\n        var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n        IPython.notebook.select(index + 1);\n    }\n}\n\nmpl.figure.prototype.handle_save = function(fig, msg) {\n    fig.ondownload(fig, null);\n}\n\n\nmpl.find_output_cell = function(html_output) {\n    // Return the cell and output element which can be found *uniquely* in the notebook.\n    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n    // IPython event is triggered only after the cells have been serialised, which for\n    // our purposes (turning an active figure into a static one), is too late.\n    var cells = IPython.notebook.get_cells();\n    var ncells = cells.length;\n    for (var i=0; i<ncells; i++) {\n        var cell = cells[i];\n        if (cell.cell_type === 'code'){\n            for (var j=0; j<cell.output_area.outputs.length; j++) {\n                var data = cell.output_area.outputs[j];\n                if (data.data) {\n                    // IPython >= 3 moved mimebundle to data attribute of output\n                    data = data.data;\n                }\n                if (data['text/html'] == html_output) {\n                    return [cell, data, j];\n                }\n            }\n        }\n    }\n}\n\n// Register the function which deals with the matplotlib target/channel.\n// The kernel may be null if the page has been refreshed.\nif (IPython.notebook.kernel != null) {\n    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n}\n"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "<div id='5f59cb94-9d3c-428d-bb7f-f605131262b3'></div>"
     },
     "metadata": {}
    },
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Progress: 98.3%"
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "interactive(children=(FloatSlider(value=0.8, description='app_use_rate', max=1.0, step=0.01), FloatSlider(valu…",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "ad244a91bf204c6ca96728e041092311"
      }
     },
     "metadata": {}
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<function __main__.update_prob(app_use_rate, report_to_app, read_notif, false_negative, time_for_test, max_notif_days, warning_after_symptoms, quarantine_after_notification)>"
     },
     "metadata": {},
     "execution_count": 11
    }
   ],
   "source": [
    "\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ((ax, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=[9,8])\n",
    "axRt = ax3.twinx()\n",
    "\n",
    "xs = []\n",
    "y_D = []\n",
    "y_MS = []\n",
    "y_MPS = []\n",
    "y_MAS = []\n",
    "y_S = []\n",
    "y_G = []\n",
    "y_Q = []\n",
    "y_InfectByASPS = []\n",
    "y_QuarantineNonI = []\n",
    "y_QuarantineI = []\n",
    "y_QuarantineNonITotal = []\n",
    "y_Test = []\n",
    "y_TestTotal = []\n",
    "y_Rt = []\n",
    "\n",
    "ax.set_ylim([0, nbIndividuals])\n",
    "\n",
    "def update_viz(graph):\n",
    "    if y_QuarantineNonITotal != []:\n",
    "        y_QuarantineNonITotal.append((graph.nbQuarantineNonI + nbIndividuals*y_QuarantineNonITotal[-1])/nbIndividuals)\n",
    "        y_TestTotal.append((graph.nbTest + nbIndividuals*y_TestTotal[-1])/nbIndividuals)\n",
    "    else:\n",
    "        y_QuarantineNonITotal.append(graph.nbQuarantineNonI/nbIndividuals)\n",
    "        y_TestTotal.append(graph.nbTest/nbIndividuals)\n",
    "\n",
    "    xs.append(len(xs))\n",
    "    y_D.append(graph.nbDead/nbIndividuals*100)\n",
    "    y_MS.append(graph.nbS/nbIndividuals*100)\n",
    "    y_MPS.append(graph.nbPS/nbIndividuals*100)\n",
    "    y_MAS.append(graph.nbAS/nbIndividuals*100)\n",
    "    y_S.append(graph.nbHealthy/nbIndividuals*100)\n",
    "    y_G.append(graph.nbCured/nbIndividuals*100)\n",
    "    y_Q.append(graph.nbQuarantineTotal)\n",
    "    y_InfectByASPS.append(graph.nbInfectedByASPS)\n",
    "    y_QuarantineNonI.append(graph.nbQuarantineNonI/nbIndividuals*100)\n",
    "    y_QuarantineI.append(graph.nbQuarantineI/nbIndividuals*100)\n",
    "    y_Test.append(graph.nbTest/nbIndividuals*100)\n",
    "\n",
    "def draw_viz(graph):\n",
    "    ax.clear()\n",
    "    ax2.clear()\n",
    "    ax3.clear()\n",
    "    ax4.clear()\n",
    "    axRt.clear()\n",
    "\n",
    "    ax.set_xlabel(\"Days\")\n",
    "    ax2.set_xlabel(\"Days\")\n",
    "    ax3.set_xlabel(\"Days\")\n",
    "    ax4.set_xlabel(\"Days\")\n",
    "\n",
    "    # computing Rt | calcul de Rt\n",
    "    for i in range(graph.stepNb):\n",
    "        if graph.contaminations[i] != 0 and graph.contaminations[i] > 5: # we just take into account days where there were more than 5 contaminations to reduce random fluctuations\n",
    "            y_Rt.append(graph.numInfectedByNewInfected[i]/graph.contaminations[i])\n",
    "        else:\n",
    "            y_Rt.append(0)\n",
    "    for i in range(1, graph.stepNb-1): # smoothing Rt curve\n",
    "        if y_Rt[i] == 0:\n",
    "            y_Rt[i] = (y_Rt[i-1] + y_Rt[i+1])/2\n",
    "\n",
    "\n",
    "    labels = [ \"Symptomatic\", \"Deceased\", \"Asymptomatic\",\"Presymptomatic\", \"Cured\", \"Healthy\"]\n",
    "    ax.stackplot(xs, y_MS, y_D, y_MAS,y_MPS, y_G, y_S, labels=labels, edgecolor=\"black\", colors=[\"red\", \"darkred\", \"orange\",\"yellow\", \"dodgerblue\", \"mediumseagreen\"])\n",
    "    ax.set_ylabel(\"Proportion of the population\")\n",
    "\n",
    "    labels2 = [\"In quarantine and non infected (percentage)\", \"In quarantine and infected (percentage)\"]\n",
    "    ax2.stackplot(xs, y_QuarantineNonI, y_QuarantineI, labels=labels2)\n",
    "    ax2.set_ylabel(\"Proportion of the population\")\n",
    "    #line, = ax3.plot(xs, y_InfectByASPS)\n",
    "    #line.set_label(\"Total infections by asympt.\")\n",
    "\n",
    "    ax3.set_ylabel(\"Quarantine days / Tests\")\n",
    "    line, = ax3.plot(xs, y_Q)\n",
    "    line.set_label(\"Cumulative quarantine days per person\")\n",
    "    line, = ax3.plot(xs, y_QuarantineNonITotal)\n",
    "    line.set_label(\"Cumulative quarantine days of healthy people per person\")\n",
    "    line, = ax3.plot(xs, y_TestTotal)\n",
    "    line.set_label(\"Cumulative number of tests per person\")\n",
    "    axRt.set_ylabel(\"Rt\", color = 'red')\n",
    "    line, = axRt.plot(xs, y_Rt, color = 'red')\n",
    "    line.set_label(\"Rt (average number of infections caused by one infected)\")\n",
    "\n",
    "    line, = ax4.plot(xs, y_Test)\n",
    "    line.set_label(\"Number of tests (in percentage of population)\")\n",
    "    ax4.set_ylabel(\"Tests\")\n",
    "\n",
    "    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=3)\n",
    "    ax2.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=1)\n",
    "    #ax3.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=1)\n",
    "    ax3.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=1)\n",
    "    #axRt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=1) #to avoid legend on top of the other\n",
    "    ax4.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), shadow=True, ncol=2)\n",
    "    plt.tight_layout()\n",
    "\n",
    "def update_prob(app_use_rate, report_to_app, read_notif, false_negative, time_for_test, max_notif_days, warning_after_symptoms, quarantine_after_notification):\n",
    "    global nbIndividuals\n",
    "    global utilApp, pReport, pReadNotif\n",
    "    global pFalseNegative, daysUntilResult, daysNotif\n",
    "    global quarantineAfterNotification, warningAfterSymptoms\n",
    "    global xs, y_D, y_MS, y_MPS, y_MAS, y_S, y_G, y_Q, y_InfectByASPS, y_Rt\n",
    "    global y_QuarantineNonI, y_QuarantineNonITotal, y_QuarantineI, y_Test, y_TestTotal\n",
    "\n",
    "    utilApp = app_use_rate\n",
    "    pReport = report_to_app\n",
    "    pReadNotif = read_notif\n",
    "    pFalseNegative = false_negative\n",
    "    daysUntilResult = time_for_test\n",
    "    daysNotif = max_notif_days\n",
    "    warningAfterSymptoms = warning_after_symptoms\n",
    "    quarantineAfterNotification = quarantine_after_notification\n",
    "    nbSteps = 60\n",
    "\n",
    "    nbIndividuals = 4000 # you may change the number of individuals for the exponential distribution graph here\n",
    "\n",
    "    sys.stdout.write('\\r' + \"Creation of the graph ...\")\n",
    "    sys.stdout.flush()\n",
    "    graph = Graph()\n",
    "    init_graph_household(graph) # default graph generation using households structure, as shown in the Results section\n",
    "    # uncomment this to get a graph with degrees following an exponential distribution\n",
    "    #init_graph_exp(graph)\n",
    "    \n",
    "    sys.stdout.write('\\r' + \"                            \")\n",
    "    sys.stdout.flush()\n",
    "    \n",
    "    xs.clear()\n",
    "    y_D.clear()\n",
    "    y_MS.clear()\n",
    "    y_MPS.clear()\n",
    "    y_MAS.clear()\n",
    "    y_S.clear()\n",
    "    y_G.clear()\n",
    "    y_Q.clear()\n",
    "    y_InfectByASPS.clear()\n",
    "    y_QuarantineNonI.clear()\n",
    "    y_QuarantineNonITotal.clear()\n",
    "    y_QuarantineI.clear()\n",
    "    y_Test.clear()\n",
    "    y_TestTotal.clear()\n",
    "    y_Rt.clear()\n",
    "\n",
    "    maxSymp = 0\n",
    "    for step_ind in range(nbSteps):\n",
    "        # update matplotlib\n",
    "        update_viz(graph)\n",
    "        # update simulation\n",
    "        step(graph)\n",
    "        sys.stdout.write('\\r'+'Progress: '+(str((100*step_ind/nbSteps))[:4]) + \"%\")\n",
    "        sys.stdout.flush()\n",
    "        maxSymp = max(maxSymp, graph.nbS)\n",
    "        \n",
    "    sys.stdout.write('\\r' + \"                     \")\n",
    "    sys.stdout.flush()\n",
    "    #usefull information to decomment if a more precise evaluation is needed\n",
    "    # print(\"Total individuals:\", nbIndividuals)\n",
    "    # print(\"Number of deceased:\", graph.nbDead)\n",
    "    # print(\"Max. nb of symptomatic people:\", maxSymp)\n",
    "    # print(\"Test per people:\", y_TestTotal[-1])\n",
    "    # print(\"Final healthy:\", y_S[-1])\n",
    "    #print(maxSymp/nbIndividuals,\",\", y_S[-1],\",\", y_Q[-1], \",\", y_TestTotal[-1])\n",
    "    draw_viz(graph)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "update_prob(utilApp, pReport, pReadNotif, pFalseNegative, daysUntilResult, daysNotif, warningAfterSymptoms, quarantineAfterNotification)\n",
    "\n",
    "interact_manual(update_prob, \\\n",
    "                 app_use_rate = widgets.FloatSlider(min=0.0, max=1.0, step=0.01, value=utilApp), \\\n",
    "                 report_to_app = widgets.FloatSlider(min=0.0, max=1.0, step=0.01, value=pReport), \\\n",
    "                 read_notif = widgets.FloatSlider(min=0.0, max=1.0, step=0.01, value=pReadNotif), \\\n",
    "                 false_negative = widgets.FloatSlider(min=0.0, max=1.0, step=0.01, value=pFalseNegative), \\\n",
    "                 time_for_test = widgets.IntSlider(min=0, max=10, step=1, value=daysUntilResult), \\\n",
    "                 max_notif_days = widgets.IntSlider(min=0, max=30, step=1, value=daysNotif), \\\n",
    "                 warning_after_symptoms = widgets.Checkbox(value=warningAfterSymptoms), \\\n",
    "                 quarantine_after_notification = widgets.Checkbox(value=quarantineAfterNotification))\n"
   ]
  },
  {
   "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.5-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}