1504 lines
63 KiB
Plaintext
Executable file
1504 lines
63 KiB
Plaintext
Executable file
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:center;\">\n",
|
|
"\n",
|
|
"<h1>\n",
|
|
"Graph project: Epidemic spreading <br/>\n",
|
|
"ENSEEIHT SN <br/>\n",
|
|
"</h1>\n",
|
|
"\n",
|
|
"<h3>\n",
|
|
"FAINSIN Laurent <br/>\n",
|
|
"GUILLOTIN Damien\n",
|
|
"</h3>"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"## Introduction\n",
|
|
" \n",
|
|
"As you know, a new epidemic has overwhelmed the world, COVID-19 jeopardizes people and changes our habits. It is easy to realise that knowing how illnesses spread is vital to our own protection. How can we predict whether a disease will cause an epidemic, how many people it will infect, which people it will infect, and whether or not it is dangerous to society as a whole ? Also, how can we determine which techniques to use in fighting an epidemic once it begins ? One way to answer all of these questions is through **mathematical modeling**.\n",
|
|
"\n",
|
|
"In this work you will have to review different epidemic modelings relying all on the representation by graphs of a human network called a **contact network**. A vertex in a contact network represents an individual and an edge between two vertices represents a contact between two individuals. The disease only spread from individual to individual if they are in contact, so through the edges. This representation is actually really common in research, and a lot of state-of-the-art modeling are built over it. From these different models you will be asked to draw conclusions from experiments on varying contact networks.\n",
|
|
"\n",
|
|
"For readability and ease of use, this project will be carried on a Jupyter Notebook, hence code and question answering have to be written in this unique file. This is a **DUO** project, no group of one person will be accepted, the duo has to be composed of same TD group students, if the number of students in the TD group is odd we will accept one group composed of three students. It will be coded in Julia using the LightGraphs package. **BEWARE:** \n",
|
|
" \n",
|
|
"- If the code does not provide good results, its readability as well as its comments are essential for the corrector to potentially find some notation points.\n",
|
|
"- The specifications of the functions have to be strictly respected.\n",
|
|
"- Do not neglect written questions they stand for an important part of the notation, you are not only evaluated on the coding. Also, even so a written question may not ask you directly to code or provide results from code, support your arguments when possible with a runable example is very welcome and sometimes even expected.\n",
|
|
"- Any initiatives and additional efforts bringing contents and thoughts out of the question scope may result in bonus points if pertinent.\n",
|
|
"\n",
|
|
"Deliverable: You will deliver your work on moodle before **23.01.2022** in a **.tar** containing the notebook with your codes and your written answers, and the different graph figures in .png you will generate. The corrector will use the student N7 computers for running your code, so take care of verifying that your work is running as expected on these computers !\n",
|
|
"\n",
|
|
"LightGraphs documentation: <https://juliagraphs.org/LightGraphs.jl/v1.2/index.html>"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"## Environment and packages installation\n",
|
|
"\n",
|
|
"**mathematical modeling**: For evaluation, coding questions have to run with no additional packages ! Only the ones present here ! However if you want to use another package to go further in your answer and add bonus contents, take care of separating the cells and precising which packages you are using."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# TO RUN ONCE\n",
|
|
"using Pkg\n",
|
|
"Pkg.activate(\".\") # Path to Manifest.toml and Project.toml\n",
|
|
"Pkg.resolve()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Import packages\n",
|
|
"using Pkg\n",
|
|
"Pkg.activate(\".\") # Path to Manifest.toml and Project.toml\n",
|
|
"\n",
|
|
"using Graphs\n",
|
|
"using Graphs: smallgraph\n",
|
|
"using Graphs: stochastic_block_model\n",
|
|
"using Graphs: barabasi_albert\n",
|
|
"using Graphs: random_regular_graph\n",
|
|
"using GraphPlot\n",
|
|
"\n",
|
|
"using Colors\n",
|
|
"using CairoMakie\n",
|
|
"using StatsBase\n",
|
|
"using Plots\n",
|
|
"using JLD2\n",
|
|
"using Compose\n",
|
|
"using Random"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"## Part 1 - SIS model\n",
|
|
"\n",
|
|
"SIS is a compartmental model, ie a model where the population is divided into subgroups that represent the disease status of its members. SIS stands for Susceptible $\\rightarrow$ Infected $\\rightarrow$ Susceptible where the susceptible group contains those who remain susceptible to the infection, and the infected group consists of those who not only have the disease but are also in the contagious period of the disease.\n",
|
|
"\n",
|
|
"\n",
|
|
"Combine with a contact network approach, this model can capture contact patterns (family, company, friends). Each vertex represents an individual in the host population, and contacts between two individuals are represented by an edge that connects the two. The probability of transmitting the disease from an infected to a susceptible individual along one of these edges or contacts is $\\beta$ (=**mathematical modeling**).\n",
|
|
"\n",
|
|
"In order for a disease to begin spreading through a network, the disease must be introduced into the population, either through infecting a proportion of the population or through infecting one individual. As time moves forward, the disease will spread away from those initially infected, and two things may occur simultaneously at each time step $t$. First, each infected individual will spread disease to each of its contacts with a probability $\\beta$. Secondly, each infectious individual will recover at a rate, $\\alpha$, at which point the individual will then no longer infect any of its contacts. After the disease has run its course, we can determine how the disease affected the network by calculating various quantities that help us better understand the outbreak.\n",
|
|
"\n",
|
|
"###### P. Van Mieghem, J. Omic, R. E. Kooij, _“Virus Spread in Networks”_, IEEE/ACM Transaction on Networking (2009)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"### 1.1 Contact networks sample"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Susceptible, Infected, Recovered, Alerted\n",
|
|
"node_colors = [colorant\"lightseagreen\", colorant\"orange\", colorant\"purple\", colorant\"red\"]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# karat7: A graph representing the karate club of N7 and the connections between the persons in this club. There are 34 people in this network. It is actually inspired by one of the most famous problem in graph theory: the Zachary's karate club.\n",
|
|
"\n",
|
|
"karat7 = smallgraph(:karate)\n",
|
|
"nodecolor = [node_colors[1]]\n",
|
|
"plot = gplot(karat7, nodefillc = nodecolor)\n",
|
|
"draw(PNG(\"img/q0_1.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# n7_2A: A graph representing the second year students at N7. Each department (SN, MF2E, 3EA) form a community where connections are denser, connections between department are rarer.\n",
|
|
"\n",
|
|
"n = [100, 70, 50]\n",
|
|
"c = [[10, 0, 0] [0.1, 10, 0] [0.1, 0.1, 10]]\n",
|
|
"n7_2A = stochastic_block_model(c, n)\n",
|
|
"nodecolor = [node_colors[1]]\n",
|
|
"plot = gplot(n7_2A, nodefillc = nodecolor)\n",
|
|
"draw(PNG(\"img/q0_2.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# toulouse_neigh: A graph representing a neighborhood composed of 1000 people in Toulouse.\n",
|
|
"\n",
|
|
"toulouse_neigh = barabasi_albert(1000, 1)\n",
|
|
"nodecolor = [node_colors[1]]\n",
|
|
"plot = gplot(toulouse_neigh, nodefillc = nodecolor)\n",
|
|
"draw(PNG(\"img/q0_3.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"### 1.2 Introduce the infection\n",
|
|
" \n",
|
|
"We denote by `state` a vector containing the disease status of each vertex where Susceptible=0 and Infected=1. Then `state` is an `Array{Int32,1}` of length the number of vertices. This array in addition of a graph (represented internally by an adjacency matrix or an adjacency list) will be the data structure of our model.\n",
|
|
" \n",
|
|
"In `Array{Int32,1}`, `Int32` refers to the kind of data in the array, here 32 bits integers, `1` refers to the dimension of the array, here we have a 1-dimensional structure so a vector"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cornflowerblue\">Question 1 (code):</strong> For each graph in the graph sample (`state`, `n7_2A`, `toulouse_neigh`) initialize the state array by assigning each vertex to susceptible and add randomly one or numerous infected people. Save the graph as a .png image using `state` and `Int32`, infected should appear in a different color (`colorant\"orange\"`).\n",
|
|
" \n",
|
|
"Due to a bug on certain version of Jupyter Notebook, the graph figures should be saved in a file and not plot inside the notebook !!!\n",
|
|
" \n",
|
|
"Gplot GitHub: <https://github.com/JuliaGraphs/GraphPlot.jl>\n",
|
|
" \n",
|
|
"Gplot examples: <https://juliagraphs.org/GraphPlot.jl/>"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"function new_state(net, ni)\n",
|
|
" return shuffle(vcat(\n",
|
|
" ones(Int, ni),\n",
|
|
" zeros(Int, nv(net) - ni)\n",
|
|
" ))\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Infect Karat7\n",
|
|
"karat7_state = new_state(karat7, 5)\n",
|
|
"nodefillc = node_colors[karat7_state.+1]\n",
|
|
"plot = gplot(karat7, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q1_1.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Infect n7_2A\n",
|
|
"n7_2A_state = new_state(n7_2A, 10)\n",
|
|
"nodefillc = node_colors[n7_2A_state.+1]\n",
|
|
"plot = gplot(n7_2A, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q1_2.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Infect Toulouse_neigh\n",
|
|
"toulouse_neigh_state = new_state(toulouse_neigh, 50)\n",
|
|
"nodefillc = node_colors[toulouse_neigh_state.+1]\n",
|
|
"plot = gplot(toulouse_neigh, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q1_3.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cornflowerblue\">Question 2 (written):</strong> What do you think/predict about the influence of the initial number of infected people and their locations on the evolution of an SIS model epidemic ?\n",
|
|
"\n",
|
|
"> Plus le nombre initial d'infectés sera élevé, et plus l'épidémie va se propager rapidement. De plus, plus la répartition initial des infectés dans la population est uniforme, plus la propagation de l'épidemie sera rapide."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"### 1.3 Spread the infection"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 3 (code):</strong> Implement the function `SIS` (respect the header and the specifications). You can use `state` to translate the probabilities. Test your algorithm on `karat7`, `n7_2A`, and `toulouse_neigh` with arbitrary $\\beta$, $\\alpha$, and $t$.\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `new_state = SIS(net,state,beta,alpha,t)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"\n",
|
|
"Take a contact network at a certain state and apply t time steps of an SIS model.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - state (Array{Int32,1}): disease status of each vertex\n",
|
|
" - beta (Float64): infection rate\n",
|
|
" - alpha (Float64): curing rate\n",
|
|
" - time (Int32): number of time step\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Int32,1}): the new state of the contact network after t time steps\n",
|
|
"\"\"\"\n",
|
|
"function SIS(net, state, beta, alpha, time)\n",
|
|
" X = zeros(time+1, 2)\n",
|
|
" adjacent_matrix = adjacency_matrix(net)\n",
|
|
" n = nv(net)\n",
|
|
"\n",
|
|
" X[1, 1] = sum(state .== 0)\n",
|
|
" X[1, 2] = sum(state .== 1)\n",
|
|
"\n",
|
|
" for t = 2:(time+1)\n",
|
|
" # println(\"-\"^50, t)\n",
|
|
" # println(\"state_0 = \", state)\n",
|
|
"\n",
|
|
" # on trouve le nombre de voisins de chaque individus\n",
|
|
" adjacent_infected = adjacent_matrix * (state .== 1)\n",
|
|
" # println(\"adjacent = \", number_infected_adjacent)\n",
|
|
"\n",
|
|
" # calcul des probas d'infection\n",
|
|
" proba_infection = (state .== 0) .* (adjacent_infected .!= 0) .* (1 .- (\n",
|
|
" (state .== 0) .* fill(1 - beta, n)\n",
|
|
" ) .^ adjacent_infected)\n",
|
|
" # println(\"proba_i = \", proba_infection\n",
|
|
"\n",
|
|
" # tirage au sort des infections\n",
|
|
" state[((state.==0).*(adjacent_infected.!=0).*rand(n)).<proba_infection] .= 1\n",
|
|
" # println(\"state_i = \", state)\n",
|
|
"\n",
|
|
" # tirage au sort des guérisons\n",
|
|
" state[((state.==1).*rand(n)).>(1-alpha)] .= 0\n",
|
|
" # println(\"state_g = \", state)\n",
|
|
"\n",
|
|
" X[t, 1] = sum(state .== 0)\n",
|
|
" X[t, 2] = sum(state .== 1)\n",
|
|
" end\n",
|
|
" return state, X\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on karat7\n",
|
|
"karat7_state = new_state(karat7, 5)\n",
|
|
"state, X = SIS(karat7, karat7_state, 0.9, 0.3, 10)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, labels=[\"Susceptible\" \"Infected\"], palette=node_colors)\n",
|
|
"png(\"img/q3_1p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(karat7, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q3_1g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on N7_2A\n",
|
|
"n7_2A_state = new_state(n7_2A, 10)\n",
|
|
"state, X = SIS(n7_2A, n7_2A_state, 0.5, 0.3, 10)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X,labels=[\"Susceptible\" \"Infected\"], palette=node_colors)\n",
|
|
"png(\"img/q3_2p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(n7_2A, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q3_2g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on Toulouse_neigh\n",
|
|
"toulouse_neigh_state = new_state(toulouse_neigh, 50)\n",
|
|
"state, X = SIS(toulouse_neigh, toulouse_neigh_state, 0.8, 0.2, 100)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, labels=[\"Susceptible\" \"Infected\"], palette=node_colors)\n",
|
|
"png(\"img/q3_3p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(toulouse_neigh, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q3_3g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"### 1.4 Simulate and understand the epidemic\n",
|
|
" \n",
|
|
"In the SIS model of this project, every disease is characterized by:\n",
|
|
"\n",
|
|
"* The infection rate $\\beta$ representing the chance of infection when being in contact with an infected individual.\n",
|
|
"* The curing rate $\\alpha$ representing the chance of being cured of the disease.\n",
|
|
"* The effective spreading rate $\\tau=\\frac{\\beta}{\\alpha}$ representing the capacity of the disease to spread. More the disease infect easily ($\\beta$ high) and less it is cured easily ($\\alpha$ low) more $\\tau$ can be high.\n",
|
|
"\n",
|
|
"We are now willing to understand what are the influences of these parameters as well as the contact network shape on an epidemic."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 4 (written):</strong> The function `SIS` you implemented launches one run of an SIS model on a given contact network. As it makes use of randomness, one run of spreading is stochastic. Then what simple method can you propose to provide a prediction of the disease spreading on a given contact network ?\n",
|
|
"\n",
|
|
"> Pour obtenir une prédiction d'un modèle SIS, on peut simplement executer la simulation un grand nombre de fois et ensuite tracer la moyenne de ces simulations."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 5 (code):</strong> Implement the function `Simulation_SIS` (respect the header and the specifications).\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `predictions, taus = Simulation_SIS(net,nbinf,betas,alphas,t,nbsimu)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"\n",
|
|
"Take a contact network, different diseases (defined by different parameters alpha and beta), a number of initial infected people and process nbsimu simulations of SIS over t time steps. You will provide the prediction of the percentage of infected at each time t as well as the spreading rate of each disease.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - nbinf (Int32): number of infected at the start of each simulation\n",
|
|
" - betas (Array{Float64,1}): array of infection rate on edges\n",
|
|
" - alphas (Array{Float64,1}): array of curing rate on vertices\n",
|
|
" - time (Int32): number of time step\n",
|
|
" - nbsimu (Int32): number of simulations\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Float64,2}): the prediction of the percentage of infected at each time step and for each disease. The first dimension contains the time steps and the second contains the diseases\n",
|
|
" - (Array{Float64,1}): effective spreading rate for each disease\n",
|
|
"\"\"\"\n",
|
|
"function Simulation_SIS(net, nbinf, betas, alphas, t, nbsimu)\n",
|
|
" taus = zeros(length(alphas))\n",
|
|
" prediction = zeros(t+1, length(taus))\n",
|
|
"\n",
|
|
" for (i, (beta, alpha)) in enumerate(zip(betas, alphas))\n",
|
|
" taus[i] = round(beta / alpha, digits = 2)\n",
|
|
" for _ in 1:nbsimu\n",
|
|
" state = new_state(net, nbinf)\n",
|
|
" _, X = SIS(net, state, beta, alpha, t)\n",
|
|
" prediction[:, i] += X[:, 2] ./ nbsimu ./ length(state)\n",
|
|
" end\n",
|
|
" end\n",
|
|
"\n",
|
|
" return prediction, taus\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 6 (written)</strong>: Run the 2 scripts below and describe what you see. Conclude on the influence of $\\tau$, $\\beta$, and $\\alpha$ on an epidemic we can model with SIS.\n",
|
|
" \n",
|
|
"> Il semblerait que plus $\\tau$ est élevé, plus le pourcentage d'infectés augmente rapidement en fonction du temps."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Script launching predictions on different diseases on karat7 and printing the precentage of infected at each time step.\n",
|
|
"betas = [0.05, 0.1, 0.01, 0.4, 0.04, 0.05, 0.005]\n",
|
|
"alphas = [0.05, 0.1, 0.01, 0.1, 0.01, 0.1, 0.01]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIS(karat7, 2, betas, alphas, 100, 1000)\n",
|
|
"plot = Plots.plot(predictions, label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q6_1.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Same as before but applied on n7_2A.\n",
|
|
"betas = [0.05, 0.1, 0.01, 0.4, 0.04, 0.05, 0.005]\n",
|
|
"alphas = [0.05, 0.1, 0.01, 0.1, 0.01, 0.1, 0.01]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIS(n7_2A, 100, betas, alphas, 100, 100)\n",
|
|
"plot = Plots.plot(predictions, label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q6_2.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Same as before but applied on toulouse_neigh. May be a bit long to run.\n",
|
|
"betas = [0.05, 0.1, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95]\n",
|
|
"alphas = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIS(toulouse_neigh, 100, betas, alphas, 100, 100)\n",
|
|
"plot = Plots.plot(predictions, label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q6_3.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 7 (written):</strong> Change the initial number of infected in the scripts above, is it in accordance with your answer in Question 2 ?\n",
|
|
"\n",
|
|
"> On remarque que plus l'on augmente le nombre initial d'infecté, plus le pourcentage d'infectés dans la population converge rapidement."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Change the initial number of infected\n",
|
|
"betas = [0.5]\n",
|
|
"alphas = [0.2]\n",
|
|
"\n",
|
|
"plot = Plots.plot()\n",
|
|
"for i in [5 10 50 100]\n",
|
|
" predictions, _ = Simulation_SIS(toulouse_neigh, i, betas, alphas, 50, 100)\n",
|
|
" Plots.plot!(predictions, label = i, xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"end\n",
|
|
"\n",
|
|
"png(\"img/q7.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cornflowerblue\">Question 8 (code):</strong> Implement a script plotting the maximum percentage of infected people according to $\\tau$ over 300 time steps for 3 contact networks:\n",
|
|
"\n",
|
|
"* A regular graph of 200 vertices with degree 2.\n",
|
|
"* A regular graph of 200 vertices with degree 5.\n",
|
|
"* A regular graph of 200 vertices with degree 10.\n",
|
|
"\n",
|
|
"You can use the function `random_regular_graph(n,d)` of LighGraphs. As you probably need to use a certain number of different values of $\\tau$ to visualize something interesting (the more there are the more the figure will be smooth) you should fix $\\alpha$ and make $\\beta$ vary. \n",
|
|
"\n",
|
|
"###### A regular graph is a graph where each vertex has the same degree."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"display(gplot(random_regular_graph(200, 2)))\n",
|
|
"display(gplot(random_regular_graph(200, 5)))\n",
|
|
"display(gplot(random_regular_graph(200, 10)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Plots of the maximum percentage of infected people according to tau over 300 time steps for 3 contact networks.\n",
|
|
"n = 200\n",
|
|
"betas = 0.01:0.01:0.99\n",
|
|
"alphas = fill(0.15, length(betas))\n",
|
|
"\n",
|
|
"plot = Plots.plot()\n",
|
|
"for d in [2 3 4 5 10]\n",
|
|
" net = random_regular_graph(n, d)\n",
|
|
" predictions, taus = Simulation_SIS(net, 5, betas, alphas, 300, 100)\n",
|
|
" max_infected = [maximum(col) for col in eachcol(predictions)]\n",
|
|
" Plots.plot!(max_infected, label = d, xlabel = \"τ\", ylabel = \"% of infected\")\n",
|
|
"end\n",
|
|
"\n",
|
|
"png(\"img/q8.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 9 (written):</strong> Describe the figure and draw conclusions on the epidemic behavior for different degrees $d$ on regular graphs. Thus, in addition of the inner properties of the disease ($\\alpha$, $\\beta$, $\\tau$) what other parameter is essential in the spreading ? Finally, what analogy can be done with real life from this experiment ?\n",
|
|
" \n",
|
|
"> Le degré $d$ du graphe représente le nombre de connections par individus. Ainsi plus $d$ est elevé plus individus se connaissent entre eux et plus l'infection est rapide. On l'observe sur le graphique, le taux d'infectés est directement proportionnel à $d$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"## Part 2 - SIR and SAIR model\n",
|
|
" \n",
|
|
"Unfortunately SIS model is valuable for diseases we can catch back since a cured person can get ill again. This is true for the flu, the cold, etc. However COVID-19 might create immunity for whom already got it and SIS can not take into account immune or dead persons. That is why we propose in this part to consider another model more adapted to COVID-19 called SIR. It stands for Susceptible $\\rightarrow$ Infected $\\rightarrow$ Recovered where the susceptible group contains those who remain susceptible to the infection, the infected group consists of those who not only have the disease but are also in the contagious period of the disease, and the recovered group contains those who were ill, got cured, are not contagious and can not get ill anymore.\n",
|
|
"\n",
|
|
"###### M. Youssef and C. Scoglio, _\"An individual-based approach to SIR epidemics in contact networks\"_, Journal of Theoretical Biology 283 (2011)\n",
|
|
" \n",
|
|
"One limitation of SIR is that it does not model the reaction of humans when they feel the presence of the epidemic. Indeed, if feeling threaten or surrounded by infected, an individual may change its behaviors: wear mask, wash its hands, etc. This result in a smaller infection rate. That is why in this part we will also consider a variant of SIR called SAIR which stands for Susceptible $\\rightarrow$ Alert $\\rightarrow$ Infected $\\rightarrow$ Recovered. A susceptible individual becomes infected by the infection rate $\\beta_0$, an infected individual recovers and gets immune by the curing rate $\\alpha$, an individual can observe the states of its neighbors, then a susceptible individual might go to the alert state if surrounded by infected individuals with an alert rate $\\kappa$ on each contact with an infected, an alert inividual becomes infected by the infection rate $\\beta_1$ where $0<\\beta_1<\\beta_0$. In our simple SAIR model, an individual can not go back to a susceptible state when he got into the alert state.\n",
|
|
"\n",
|
|
"###### F. Darabi Sahneh and C. Scoglio, _\"Epidemic Spread in Human Networks\"_, 50th IEEE Conf. Decision and Contol, Orlando, Florida (2011)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"### 2.1 SIR\n",
|
|
" \n",
|
|
"The vector containing the disease status `state` has to change a bit since we added a new state. Hence it will be an `Array{Int32,1}` where Susceptible=0, Infected=1, and Recovered=2."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 10 (code):</strong> Implement the function `SIR` (respect the header and the specifications). You can use `state` to translate the probabilities. Test your algorithm on `state`, `n7_2A`, and `toulouse_neigh` with arbitrary $\\beta$, $\\alpha$, and $t$. Recovered vertices should appear in a different color (`colorant\"purple\"`).\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `new_state = SIR(net,state,beta,alpha,t)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"\n",
|
|
"Take a contact network at a certain state and apply t time steps of an SIR model.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - state (Array{Int32,1}): disease status of each vertex\n",
|
|
" - beta (Float64): infection rate\n",
|
|
" - alpha (Float64): curing rate\n",
|
|
" - time (Int32): number of time step\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Int32,1}): The new state of the contact network after t time steps\n",
|
|
"\"\"\"\n",
|
|
"function SIR(net, state, beta, alpha, time)\n",
|
|
" X = zeros(time+1, 3)\n",
|
|
" adjacent_matrix = adjacency_matrix(net)\n",
|
|
" n = nv(net)\n",
|
|
"\n",
|
|
" X[1, 1] = sum(state .== 0)\n",
|
|
" X[1, 2] = sum(state .== 1)\n",
|
|
" X[1, 3] = sum(state .== 2)\n",
|
|
"\n",
|
|
" for t = 2:(time+1)\n",
|
|
" # println(\"-\"^50, t)\n",
|
|
" # println(\"state_0 = \", state)\n",
|
|
"\n",
|
|
" # on trouve le nombre de voisins de chaque individus\n",
|
|
" adjacent_infected = adjacent_matrix * (state .== 1)\n",
|
|
" # println(\"adjacent = \", number_infected_adjacent)\n",
|
|
"\n",
|
|
" # calcul des probas d'infection\n",
|
|
" proba_infection = (state .== 0) .* (adjacent_infected .!= 0) .* (1 .- (\n",
|
|
" (state .== 0) .* fill(1 - beta, n)\n",
|
|
" ) .^ adjacent_infected)\n",
|
|
" # println(\"proba_i = \", proba_infection\n",
|
|
"\n",
|
|
" # tirage au sort des infections\n",
|
|
" state[((state.==0).*(adjacent_infected.!=0).*rand(n)).<proba_infection] .= 1\n",
|
|
" # println(\"state_i = \", state)\n",
|
|
"\n",
|
|
" # tirage au sort des guérisons\n",
|
|
" state[((state.==1).*rand(n)).>=(1-alpha)] .= 2\n",
|
|
" # println(\"state_g = \", state)\n",
|
|
"\n",
|
|
" X[t, 1] = sum(state .== 0)\n",
|
|
" X[t, 2] = sum(state .== 1)\n",
|
|
" X[t, 3] = sum(state .== 2)\n",
|
|
" end\n",
|
|
" return state, X\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on Karat7\n",
|
|
"karat7_state = new_state(karat7, 1)\n",
|
|
"state, X = SIR(karat7, karat7_state, 0.7, 0.1, 10)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, labels=[\"Susceptible\" \"Infected\" \"Recovered\"], palette=node_colors)\n",
|
|
"png(\"img/q10_1p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(karat7, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q10_1g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on n7_2A\n",
|
|
"n7_2A_state = new_state(n7_2A, 5)\n",
|
|
"state, X = SIR(n7_2A, n7_2A_state, 0.9, 0.1, 40)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, labels=[\"Susceptible\" \"Infected\" \"Recovered\"], palette=node_colors)\n",
|
|
"png(\"img/q10_2p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(n7_2A, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q10_2g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on toulouse_neigh\n",
|
|
"toulouse_neigh_state = new_state(toulouse_neigh, 50)\n",
|
|
"state, X = SIR(toulouse_neigh, toulouse_neigh_state, 0.9, 0.1, 40)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, labels=[\"Susceptible\" \"Infected\" \"Recovered\"], palette=node_colors)\n",
|
|
"png(\"img/q10_3p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(toulouse_neigh, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q10_3g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 11 (code):</strong> Implement the function `Simulation_SIR` (respect the header and the specifications).\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `predictions, taus = Simulation_SIR(net,nbinf,betas,alphas,t,nbsimu)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"Take a contact network, different diseases (defined by different parameters alpha and beta), a number of initial infected people and process nbsimu simulations of SIR over t time steps. You will provide the prediction of the percentage of infected at each time t as well as the spreading rate of each disease.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - nbinf (Int32): number of infected at the start of each simulation\n",
|
|
" - betas (Array{Float64,1}): array of infection rate on edges\n",
|
|
" - alphas (Array{Float64,1}): array of curing rate on vertices\n",
|
|
" - t (Int32): number of time step\n",
|
|
" - nbsimu (Int32): number of simulations\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Float64,3}): the prediction of the percentage of infected, the percentage of susceptible and the percentage of recovered at each time step and for each disease. The first dimension contains the time steps,the second contains the diseases, and the third the status(Infected: [:,:,1], Recovered: [:,:,2], Susceptible: [:,:,3])\n",
|
|
" - (Array{Float64,1}): effective spreading rate for each disease\n",
|
|
"\"\"\"\n",
|
|
"function Simulation_SIR(net, nbinf, betas, alphas, t, nbsimu)\n",
|
|
" taus = zeros(length(alphas))\n",
|
|
" prediction = zeros(t+1, length(taus), 3)\n",
|
|
"\n",
|
|
" for (i, (beta, alpha)) in enumerate(zip(betas, alphas))\n",
|
|
" taus[i] = round(beta / alpha, digits = 2)\n",
|
|
" for _ in 1:nbsimu\n",
|
|
" state = new_state(net, nbinf)\n",
|
|
" _, X = SIR(net, state, beta, alpha, t)\n",
|
|
" prediction[:, i, :] += X ./ nbsimu ./ length(state)\n",
|
|
" end\n",
|
|
" end\n",
|
|
"\n",
|
|
" return prediction, taus\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 12 (written):</strong> Run the script below and describe what you see. Why the infected curve does not behave the same as for SIS ? \n",
|
|
" \n",
|
|
"> TODO"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Script launching prediction on one disease on n7_2A and plotting the percentage\n",
|
|
"# of infected, susceptible and recovered at each time step.\n",
|
|
"predictions, taus = Simulation_SIR(n7_2A, 2, [0.3], [0.2], 50, 1000)\n",
|
|
"\n",
|
|
"plot = Plots.plot(\n",
|
|
" [predictions[:, :, 1] predictions[:, :, 2] predictions[:, :, 3]],\n",
|
|
" label = [\"Infected\" \"Recovered\" \"Susceptible\"],\n",
|
|
" xlabel = \"t\",\n",
|
|
" ylabel = \"%\",\n",
|
|
" palette=node_colors\n",
|
|
")\n",
|
|
"png(\"img/q12.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 13 (written):</strong> As for Question 6 script 2 plot the evolution of the percentage of infected for many $\\tau$. Describe what you see.\n",
|
|
"\n",
|
|
"> TODO"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Equivalent experiment as for Question 6 script 2\n",
|
|
"# Script launching predictions on different diseases on karat7 and printing the precentage of infected at each time step.\n",
|
|
"betas = [0.05, 0.1, 0.01, 0.4, 0.04, 0.05, 0.005]\n",
|
|
"alphas = [0.05, 0.1, 0.01, 0.1, 0.01, 0.1, 0.01]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIR(karat7, 2, betas, alphas, 100, 1000)\n",
|
|
"plot = Plots.plot(predictions[:, :, 2], label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q13_1.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"betas = [0.05, 0.1, 0.01, 0.4, 0.04, 0.05, 0.005]\n",
|
|
"alphas = [0.05, 0.1, 0.01, 0.1, 0.01, 0.1, 0.01]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIR(n7_2A, 10, betas, alphas, 100, 1000)\n",
|
|
"plot = Plots.plot(predictions[:, :, 2], label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q13_2.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Same as before but applied on toulouse_neigh. May be a bit long to run.\n",
|
|
"betas = [0.05, 0.1, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95]\n",
|
|
"alphas = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]\n",
|
|
"\n",
|
|
"predictions, taus = Simulation_SIR(toulouse_neigh, 100, betas, alphas, 100, 1000)\n",
|
|
"plot = Plots.plot(predictions[:, :, 2], label = taus', xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"png(\"img/q13_3.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cornflowerblue\">Question 14 (code):</strong> Implement a script plotting the number of infected over 75 time steps for $\\beta=0.3$ and $\\alpha=0.2$ fixed and on 3 contact networks:\n",
|
|
" \n",
|
|
"* A regular graph of 200 vertices with degree 2.\n",
|
|
"* A regular graph of 200 vertices with degree 5.\n",
|
|
"* A regular graph of 200 vertices with degree 10."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Plots of the number of infected people according to tau over 75 time steps for 3 contact networks.\n",
|
|
"n = 200\n",
|
|
"Plots.plot()\n",
|
|
"betas = [0.3]\n",
|
|
"alphas = [0.2]\n",
|
|
"\n",
|
|
"plot = Plots.plot()\n",
|
|
"for d in [2 5 10]\n",
|
|
" net = random_regular_graph(n, d)\n",
|
|
" predictions, taus = Simulation_SIR(net, 1, betas, alphas, 75, 1000)\n",
|
|
" Plots.plot!(predictions[:, :, 2], label = d, xlabel = \"Time\", ylabel = \"% of infected\")\n",
|
|
"end\n",
|
|
"Plots.plot!()\n",
|
|
"\n",
|
|
"png(\"img/q14.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 15 (written):</strong> From the previous figure, explain why lockdown can be interesting when hospital places are lacking ?\n",
|
|
"\n",
|
|
"> On observe qu'en diminuant le nombre de contacts potentiels entre les personnes d'un réseau, le nombre maximal d'infecté à un temps donné diminue. Ainsi par la mise en place d'un confinement les hopitaux peuvent être capable de supporter plus facilement la charge d'infectés."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"### 2.2 SAIR\n",
|
|
" \n",
|
|
"The vector containing the disease status `state` has to change a bit since we added a new state. Hence it will be an `Array{Int32,1}` where Susceptible=0, Infected=1, Recovered=2, and Alert=3."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 16 (code):</strong> Implement the function `SAIR` (respect the header and the specifications). You can use `state` to translate the probabilities. Test your algorithm on `state`, `n7_2A`, and `toulouse_neigh` with arbitrary $\\beta$, $\\alpha$, and $t$. Alerted vertices should appear in a different color (`colorant\"lightgreen\"`).\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `new_state = SAIR(net,state,beta0,beta1,alpha,kappa,t)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"Take a contact network at a certain state and apply t time steps of an SAIR model.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - state (Array{Int32,1}): disease status of each vertex\n",
|
|
" - beta0 (Float64): infection rate when not alert\n",
|
|
" - beta1 (Float64): infection rate when alert\n",
|
|
" - alpha (Float64): curing rate\n",
|
|
" - kappa (Float64): alerting rate\n",
|
|
" - t (Int32): number of time step\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Int32,1}): The new state of the contact network after t time steps\n",
|
|
"\"\"\"\n",
|
|
"function SAIR(net, state, beta0, beta1, alpha, kappa, time)\n",
|
|
" X = zeros(time+1, 4)\n",
|
|
" adjacent_matrix = adjacency_matrix(net)\n",
|
|
" n = nv(net)\n",
|
|
"\n",
|
|
" X[1, 1] = sum(state .== 0)\n",
|
|
" X[1, 2] = sum(state .== 1)\n",
|
|
" X[1, 3] = sum(state .== 2)\n",
|
|
" X[1, 4] = sum(state .== 3)\n",
|
|
"\n",
|
|
" for t = 2:(time+1)\n",
|
|
" # println(\"-\"^50, t)\n",
|
|
" # println(\"state_0 = \", state)\n",
|
|
"\n",
|
|
" # on trouve le nombre de voisins de chaque individus\n",
|
|
" adjacent_infected = adjacent_matrix * (state .== 1)\n",
|
|
" # println(\"adjacent = \", number_infected_adjacent)\n",
|
|
"\n",
|
|
" # calcul des probas d'alert\n",
|
|
" proba_alert = (state .== 0) .* (adjacent_infected .!= 0) .* (1 .- (\n",
|
|
" (state .== 0) .* fill(1 - kappa, n)\n",
|
|
" ) .^ adjacent_infected)\n",
|
|
" # println(\"proba_a = \", proba_alert)\n",
|
|
"\n",
|
|
" # tirage au sort les alerts\n",
|
|
" state[((state.==0).*(adjacent_infected.!=0).*rand(n)).<proba_alert] .= 3\n",
|
|
" # println(\"state_a = \", state)\n",
|
|
"\n",
|
|
" # calcul des probas d'infection\n",
|
|
" proba_infection = ((state .== 0) .+ (state .== 3)) .* (adjacent_infected .!= 0) .* (1 .- (\n",
|
|
" (state .== 0) .* fill(1 - beta0, n)\n",
|
|
" .+\n",
|
|
" (state .== 3) .* fill(1 - beta1, n)\n",
|
|
" ) .^ adjacent_infected)\n",
|
|
" # println(\"proba_i = \", proba_infection)\n",
|
|
"\n",
|
|
" # tirage au sort des infections\n",
|
|
" state[(((state.==0).+(state.==3)).*(adjacent_infected.!=0).*rand(n)).<proba_infection] .= 1\n",
|
|
" # println(\"state_i = \", state)\n",
|
|
"\n",
|
|
" # tirage au sort des guérisons\n",
|
|
" state[((state.==1).*rand(n)).>=(1-alpha)] .= 2\n",
|
|
" # println(\"state_g = \", state)\n",
|
|
"\n",
|
|
" X[t, 1] = sum(state .== 0)\n",
|
|
" X[t, 2] = sum(state .== 1)\n",
|
|
" X[t, 3] = sum(state .== 2)\n",
|
|
" X[t, 4] = sum(state .== 3)\n",
|
|
" end\n",
|
|
" return state, X\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on Karat7\n",
|
|
"karat7_state = new_state(karat7, 5)\n",
|
|
"state, X = SAIR(karat7, karat7_state, 0.9, 0.1, 0.1, 0.5, 10)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, label = [\"Susceptible\" \"Infected\" \"Recovered\" \"Alert\"], palette=node_colors)\n",
|
|
"png(\"img/q16_1p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(karat7, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q16_1g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on N7_2A\n",
|
|
"n7_2A_state = new_state(n7_2A, 10)\n",
|
|
"state, X = SAIR(n7_2A, n7_2A_state, 0.9, 0.05, 0.1, 0.5, 30)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, label = [\"Infected\" \"Recovered\" \"Susceptible\" \"Alert\"], palette=node_colors)\n",
|
|
"png(\"img/q16_2p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(n7_2A, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q16_2g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Test on Toulouse_neigh\n",
|
|
"toulouse_neigh_state = new_state(toulouse_neigh, 50)\n",
|
|
"state, X = SAIR(toulouse_neigh, toulouse_neigh_state, 0.9, 0.1, 0.1, 0.5, 50)\n",
|
|
"\n",
|
|
"plot = Plots.plot(1:size(X)[1], X, label = [\"Infected\" \"Recovered\" \"Susceptible\" \"Alert\"], palette=node_colors)\n",
|
|
"png(\"img/q16_3p.png\")\n",
|
|
"display(plot)\n",
|
|
"\n",
|
|
"nodefillc = node_colors[state.+1]\n",
|
|
"plot = gplot(toulouse_neigh, nodefillc = nodefillc)\n",
|
|
"draw(PNG(\"img/q16_3g.png\", 100cm, 100cm), plot)\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 17 (code):</strong> Implement the function `Simulation_SAIR` (respect the header and the specifications).\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `predictions, taus = Simulation_SAIR(net,nbinf,betas0,betas1,alphas,kappas,t,nbsimu)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"\n",
|
|
"Take a contact network, different diseases (defined by different parameters alpha and beta), a number of initial infected people and process nbsimu simulations of SAIR over t time steps. You will provide the prediction of the percentage of infected at each time t as well as the spreading rate of each disease.\n",
|
|
"\n",
|
|
"* PARAMS\n",
|
|
" - net (LightGraph): graph representing the contact network\n",
|
|
" - nbinf (Int32): number of infected at the start of each simulation\n",
|
|
" - betas0 (Array{Float64,1}): array of infection rate when not alert on edges\n",
|
|
" - betas1 (Array{Float64,1}): array of infection rate when alert on edges\n",
|
|
" - alphas (Array{Float64,1}): array of curing rate on vertices\n",
|
|
" - kappas (Array{Float64,1}): array of alerting rate on edges\n",
|
|
" - t (Int32): number of time step\n",
|
|
" - nbsimu (Int32): number of simulations\n",
|
|
"\n",
|
|
"* RETURNS\n",
|
|
" - (Array{Float64,3}): the prediction of the percentage of infected, the percentage of susceptible and the percentage of recovered at each time step and for each disease. The first dimension contains the time steps, the second contains the diseases, and the third the status (Infected: [:,:,1], Recovered: [:,:,2], Susceptible: [:,:,3])\n",
|
|
" - (Array{Float64,1}): effective spreading rate for each disease\n",
|
|
"\"\"\"\n",
|
|
"function Simulation_SAIR(net, nbinf, betas0, betas1, alphas, kappas, t, nbsimu)\n",
|
|
" taus = zeros(length(alphas))\n",
|
|
" prediction = zeros(t+1, length(taus), 4)\n",
|
|
"\n",
|
|
" for (i, (beta0, beta1, alpha, kappa)) in enumerate(zip(betas0, betas1, alphas, kappas))\n",
|
|
" taus[i] = round(beta0 / alpha, digits = 2)\n",
|
|
" for _ in 1:nbsimu\n",
|
|
" state = new_state(net, nbinf)\n",
|
|
" _, X = SAIR(net, state, beta0, beta1, alpha, kappa, t)\n",
|
|
" prediction[:, i, :] += X ./ nbsimu ./ length(state)\n",
|
|
" end\n",
|
|
" end\n",
|
|
"\n",
|
|
" return prediction, taus\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 18 (written):</strong> Run the script below comparing the number of infected of SAIR and SIR and comment what you see.\n",
|
|
"\n",
|
|
"> TODO"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Script launching prediction on one disease on toulouse_neigh and plotting the percentage\n",
|
|
"# of infected at each time step for SIR and SAIR.\n",
|
|
"predictions1, _ = Simulation_SAIR(toulouse_neigh, 2, [0.2], [0.1], [0.1], [0.5], 100, 1000)\n",
|
|
"predictions2, _ = Simulation_SIR( toulouse_neigh, 2, [0.2], [0.1], 100, 1000)\n",
|
|
"\n",
|
|
"plot = Plots.plot(\n",
|
|
" [predictions1[:, :, 2] predictions2[:, :, 2]],\n",
|
|
" label = [\"SAIR\" \"SIR\"],\n",
|
|
" xlabel = \"t\",\n",
|
|
" ylabel = \"%\"\n",
|
|
")\n",
|
|
"png(\"img/q18.png\")\n",
|
|
"display(plot)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 19 (written):</strong> Of course the presented SIS, SIR, and SAIR models are limitated in their modelization of the reality. Formulate few of these limitations (at least 2). Propose few algorithm addons/ideas (at least 2) which would make the models more complex and more accurate in regards to the reality.\n",
|
|
" \n",
|
|
"> TODO"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"## Part 3 - Discover patient zero\n",
|
|
" \n",
|
|
"In the two previous parts you may have realised that understanding and controlling the spread of epidemics on contact networks is an important task. However, information about the origin of the epidemic could be also extremely useful to reduce or prevent future outbreaks. Thus, in this part we will focus on algorithm solutions to answer this issue.\n",
|
|
" \n",
|
|
"The stochastic nature of infection propagation makes the estimation of the epidemic origin intrinsically hard: indeed, different initial conditions can lead to the same configuration at the observation time. Methods such as the distance centrality or the Jordan center try to approximate it. They both rely on spatial information by stating that the first infected is probably at the center of the cluster of infection. Mathematically:\n",
|
|
" \n",
|
|
"* The jordan center is expressed as $\\min_{v\\in \\mathcal{I}}\\max_{n\\in \\mathcal{I}}d(v,n)$ where $\\mathcal{I}$ is a connected component of the original contact network containing all infected and recovered vertices, and where $d(\\cdot,\\cdot)$ is the distance (= the shortest path) between 2 vertices (if not weighted graph each edge accounts for 1 unit). \n",
|
|
"* The distance centrality is expressed as $\\min_{v\\in \\mathcal{I}}\\sum_{n\\in \\mathcal{I}}d(v,n)(\\delta_{n,I} + \\delta_{n,R}/\\alpha)$, where $\\delta_{n,I}=1$ if the vertex n is infected ($=0$ otherwise), and where $\\delta_{n,R}=1$ if the vertex n is recovered ($=0$ otherwise). You may note that in distance centrality we increase the weight of the recovered vertices by a factor $1/\\alpha$, it translates the fact that recovered vertices tend to be closer to the origin of the epidemic since they probably got ill before.\n",
|
|
" \n",
|
|
" \n",
|
|
"We formulate the problem as follow: given a contact network and a snapshot of epidemic spread at a certain time, determine the infection source. A snapshot is a given `state` array for a contact network.\n",
|
|
"\n",
|
|
"###### A. Y. Lokhov, M. Mézard, H. Ohta, and L. Zdeborová, _\"Inferring the origin of an epidemic with a dynamic message-passing algorithm\"_, Physical Review (2014)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 20 (code):</strong> Implement the function `jordan` (respect the header and the specifications). You will need to use the function `dijkstra_shortest_paths` of the LighGraphs library, refer to the doc for more information. If there are multiple minimal vertices, then return the first one.\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `zero = jordan(g,state,alpha)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"function jordan(g, state)\n",
|
|
" \"\"\"Find patient zero by mean of the jordan center method.\n",
|
|
"\n",
|
|
" PARAMS\n",
|
|
" g (LightGraph): graph representing the contact network\n",
|
|
" state (Array{Int32,1}): disease status of each vertex\n",
|
|
"\n",
|
|
" RETURNS\n",
|
|
" (Int32): the patient zero vertex number \n",
|
|
" \"\"\"\n",
|
|
" # TODO\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
" \n",
|
|
"<strong style=\"color:cornflowerblue\">Question 21 (code):</strong> Implement the function `distance` (respect the header and the specifications). You will need to use the function `dijkstra_shortest_paths` of the LighGraphs library, refer to the doc for more information. If there are multiple minimal vertices, then return the first one.\n",
|
|
" \n",
|
|
"###### The corrector should be able to write `zero = distance(g,state,alpha)` with your code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"function distance(g, state, alpha = 1.0)\n",
|
|
" \"\"\"Find patient zero by mean of the distance centrality method.\n",
|
|
"\n",
|
|
" PARAMS\n",
|
|
" g (LightGraph): graph representing the contact network\n",
|
|
" state (Array{Int32,1}): disease status of each vertex\n",
|
|
" alpha (Float64): curing rate\n",
|
|
"\n",
|
|
" RETURNS\n",
|
|
" (Int32): the patient zero vertex number \n",
|
|
" \"\"\"\n",
|
|
" # TODO\n",
|
|
"end"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<div style=\"width:75%;margin:0 auto;text-align:justify;\">\n",
|
|
"\n",
|
|
"<strong style=\"color:cadetblue\">Question 22 (written):</strong> Run the 3 following scripts using your functions `state` and `state` and comment on the results.\n",
|
|
" \n",
|
|
"The contact network is karat7 for 2 different patient zero and a $50\\times 50$ grid. The real patient zero (\"Z\"), your jordan (\"J\") and distance (\"D\") approximations are appearing in `colorant\"lightblue\"`.\n",
|
|
"\n",
|
|
"> TODO"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Loading a snapshot of karat7\n",
|
|
"@load \"karat7_Q22_1.jld2\" g state pat_zero alpha beta loc_x loc_y\n",
|
|
"\n",
|
|
"# Run the patient zero finding function\n",
|
|
"cent_pat_zero = distance(g, state, alpha)\n",
|
|
"jor_pat_zero = jordan(g, state)\n",
|
|
"\n",
|
|
"# Some display options \n",
|
|
"labels = Array{String,1}(undef, nv(g))\n",
|
|
"for k = 1:nv(g)\n",
|
|
" if state[k] == 1\n",
|
|
" labels[k] = \"I\"\n",
|
|
" elseif state[k] == 2\n",
|
|
" labels[k] = \"R\"\n",
|
|
" else\n",
|
|
" labels[k] = \"S\"\n",
|
|
" end\n",
|
|
"end\n",
|
|
"\n",
|
|
"if cent_pat_zero == jor_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J+Z\"\n",
|
|
"elseif cent_pat_zero == jor_pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"elseif cent_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+Z\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
"elseif jor_pat_zero == pat_zero\n",
|
|
" labels[jor_pat_zero] = \"J+Z\"\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
"else\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"end\n",
|
|
"\n",
|
|
"nodecolor = [colorant\"lightseagreen\", colorant\"orange\", colorant\"purple\"]\n",
|
|
"colors = nodecolor[state+ones(Int32, nv(g))]\n",
|
|
"colors[pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[cent_pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[jor_pat_zero] = colorant\"lightblue\"\n",
|
|
"\n",
|
|
"# Display\n",
|
|
"draw(PNG(\"karat7_Q22_1.png\", 20cm, 20cm), gplot(g, loc_x, loc_y, nodefillc = colors, nodelabel = labels))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Loading a snapshot of karat7\n",
|
|
"@load \"karat7_Q22_2.jld2\" g state pat_zero alpha beta loc_x loc_y\n",
|
|
"\n",
|
|
"# Run the patient zero finding function\n",
|
|
"cent_pat_zero = distance(g, state, alpha)\n",
|
|
"jor_pat_zero = jordan(g, state)\n",
|
|
"\n",
|
|
"# Some display options \n",
|
|
"labels = Array{String,1}(undef, nv(g))\n",
|
|
"for k = 1:nv(g)\n",
|
|
" if state[k] == 1\n",
|
|
" labels[k] = \"I\"\n",
|
|
" elseif state[k] == 2\n",
|
|
" labels[k] = \"R\"\n",
|
|
" else\n",
|
|
" labels[k] = \"S\"\n",
|
|
" end\n",
|
|
"end\n",
|
|
"\n",
|
|
"if cent_pat_zero == jor_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J+Z\"\n",
|
|
"elseif cent_pat_zero == jor_pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"elseif cent_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+Z\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
"elseif jor_pat_zero == pat_zero\n",
|
|
" labels[jor_pat_zero] = \"J+Z\"\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
"else\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"end\n",
|
|
"\n",
|
|
"nodecolor = [colorant\"lightseagreen\", colorant\"orange\", colorant\"purple\"]\n",
|
|
"colors = nodecolor[state+ones(Int32, nv(g))]\n",
|
|
"colors[pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[cent_pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[jor_pat_zero] = colorant\"lightblue\"\n",
|
|
"\n",
|
|
"# Display\n",
|
|
"draw(PNG(\"karat7_Q22_2.png\", 20cm, 20cm), gplot(g, loc_x, loc_y, nodefillc = colors, nodelabel = labels))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Loading a snapshot of grid50\n",
|
|
"@load \"grid50_Q22.jld2\" g state pat_zero alpha beta loc_x loc_y\n",
|
|
"\n",
|
|
"# Run the patient zero finding function\n",
|
|
"cent_pat_zero = distance(g, state, alpha)\n",
|
|
"jor_pat_zero = jordan(g, state)\n",
|
|
"\n",
|
|
"# Some display options \n",
|
|
"labels = Array{String,1}(undef, nv(g))\n",
|
|
"for k = 1:nv(g)\n",
|
|
" if state[k] == 1\n",
|
|
" labels[k] = \"I\"\n",
|
|
" elseif state[k] == 2\n",
|
|
" labels[k] = \"R\"\n",
|
|
" else\n",
|
|
" labels[k] = \"S\"\n",
|
|
" end\n",
|
|
"end\n",
|
|
"\n",
|
|
"if cent_pat_zero == jor_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J+Z\"\n",
|
|
"elseif cent_pat_zero == jor_pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"elseif cent_pat_zero == pat_zero\n",
|
|
" labels[cent_pat_zero] = \"C+Z\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
"elseif jor_pat_zero == pat_zero\n",
|
|
" labels[jor_pat_zero] = \"J+Z\"\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
"else\n",
|
|
" labels[cent_pat_zero] = \"C\"\n",
|
|
" labels[jor_pat_zero] = \"J\"\n",
|
|
" labels[pat_zero] = \"Z\"\n",
|
|
"end\n",
|
|
"\n",
|
|
"nodecolor = [colorant\"lightseagreen\", colorant\"orange\", colorant\"purple\"]\n",
|
|
"colors = nodecolor[state+ones(Int32, nv(g))]\n",
|
|
"colors[pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[cent_pat_zero] = colorant\"lightblue\"\n",
|
|
"colors[jor_pat_zero] = colorant\"lightblue\"\n",
|
|
"\n",
|
|
"# Display\n",
|
|
"draw(PNG(\"grid50_Q22.png\", 100cm, 100cm), gplot(g, loc_x, loc_y, nodefillc = colors, nodelabel = labels))"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Julia 1.8.0-DEV",
|
|
"language": "julia",
|
|
"name": "julia-1.8"
|
|
},
|
|
"language_info": {
|
|
"file_extension": ".jl",
|
|
"mimetype": "application/julia",
|
|
"name": "julia",
|
|
"version": "1.8.0-DEV"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|