用Julia模拟图形上的简单流行性传播过程
图是表示一组对象和它们之间关系的结构。
从形式上看，我们说对象是顶点或节点，关系是边或链接。我们经常可以看到图出现在我们周围。
比如说。
- 顶点是城市，边是直接连接两个城市的道路。
- 顶点是网络上的网页，边是连接两个网页的链接。
- 顶点是人，边是两个人之间的友谊。
这些只是图的一些应用--这个清单非常庞大。
我们通常用视觉来表示图，把顶点画成圆，把边画成连接两个圆的线。
举例来说。
这是一个顶点为{1, 2, 3, 4, 5, 6} ，边为{(1, 2), (1, 5), (2, 5), (2, 3), (3, 4), (4, 5), (4, 6)} 的图形。
在本教程中，我们将使用图来模拟流行病的传播。
前提条件
对于本教程，读者将需要。
- 安装了Julia。
- 安装了Jupyter-Lab。
- 安装FFmpeg并添加到Path中（可选--仅用于生成动画）。
- 基本的编程知识。
- 图形的概念。
如果你也有一些Julia的背景就更好了，但如果你没有，那也不是问题如果你有一些编程知识，你就不会有很大的困难去理解这些代码片段在做什么。
扩散的动态
由于图的性质，即表示成对的相互作用/关系，它们是模拟流行病传播的良好结构。我们可以把个体（例如人或动物）表示为顶点，把接触/接触情况表示为边。
我们将对该流行病进行如下建模。
- 每个顶点都将处于以下状态之一。
- S - 易感者
- I - 感染者
- R - 康复的
- 在每一个时间步骤中。
- 一些暴露的易感顶点（以前被感染顶点的邻居）将被感染。
- 之前被感染的顶点会恢复，不再是易感或感染者。
这是个简单的模型，实际上也是为了这个目的!但在本教程结束时，你将能够根据自己的需要进行修改并添加新的功能或动态。
第1步：设置好东西
首先，让我们来安装我们需要的模块。对于这个项目，我们将需要以下包。
打开你的终端并输入。
> julia
> using Pkg
> Pkg.add("Graphs")
> Pkg.add("MetaGraphs")
> Pkg.add("Plots")
> Pkg.add("GraphPlot")
> Pkg.add("StatsBase")
> Pkg.add("Cairo")
> Pkg.add("Fontconfig")
> Pkg.add("Compose")
> exit()
安装完软件包后，再打开Jupyter-lab。
> jupyter-lab
为本教程创建一个新的Julia笔记本，让我们进入下一个步骤
第2步：创建并绘制一个简单的图形
如果我们只需要操作简单的图形，Graphs.jl包就足够了。但是，正如我们前面看到的，每个顶点都需要一些元数据：它的状态。所以我们还需要MetaGraphs.jl包。让我们来导入它们。
using Graphs, MetaGraphs
创建你的第一个图。
G = MetaGraph()
你会看到这个图是一个{0, 0} undirected Int64 metagraph 。这意味着G 有0 的顶点和0 的边，它是无方向的--它的边不是*"单向的"。*现在让我们添加一些顶点和一些边。
# adding 6 new vertices
add_vertices!(G, 6)
# creating edges
add_edge!(G, 1, 2)
add_edge!(G, 1, 3)
add_edge!(G, 2, 3)
add_edge!(G, 3, 4)
add_edge!(G, 4, 5)
add_edge!(G, 5, 6)
注意在函数名称前面的
!。这是个惯例。我们说带有!的函数会修改其参数。在这种情况下，add_vertices!和add_edge!通过添加新的顶点和新的边来修改图形G。
让我们来看看我们的图。导入GraphPlot 。
using GraphPlot
绘制图形是非常简单的。
只要运行。
# plotting the Graph
gplot(G)
你应该看到像这样的东西*（也许不完全是）*。
我们的第一个图表!
第3步：设置和获取道具
到目前为止，我们只有一个简单的图形，没有其他信息。v 我们通过向set_props! 传递一个包含所需信息的字典来向顶点添加元数据。让我们写一个函数，将G 中的所有顶点的状态设置为S （易受影响）。
function init_nodes_states!(G)
for v in vertices(G)
set_props!(G, v, Dict(:state=>"S", :color=>"yellow"))
end
end
通过init_nodes_states! ，我们将图G 的所有顶点的状态设置为S （易感）。
请注意，我们也在添加信息
color。那是为了可视化的目的。让我们这样约定。
- 黄色：易受感染
- 红色：被感染
- 蓝色：已恢复
让我们把init_nodes_states! 到G 。
init_nodes_states!(G)
要检查一个给定顶点的所有道具，比方说顶点1，我们可以运行。
# checking a prop for a vertex
MetaGraphs.props(G, 1)
这应该会返回。
Dict{Symbol,Any} with 2 entries:
:color => "yellow"
:state => "S"
我们还可以使用get_prop ，获得一个特定的道具，比如说color 。
get_prop(G, 1, :color)
这段代码从图形G ，返回顶点1 的道具:color 。
让我们用新的颜色绘制我们的图。
function plot_graph(G)
fillcolors = [get_prop(G, i, :color) for i in vertices(G)]
return gplot(G, nodefillc=fillcolors, layout=spectral_layout)
end
plot_graph(G)
在plot_graph 中，我们正在为G 中的每个顶点计算color 道具，并将颜色传递给gplot 的参数nodefillc 。我们现在也在使用spectral_layout 。根据NetworkLayout.jl网站，这是*"一种未被重视的图形布局方法；比更常见的基于弹簧的方法更容易、更简单、更快速"。此外，标准的（基于弹簧的）*布局方法每次都会生成不同的图形，这对我们的可视化效果不好。
第4步：动态编程
我们现在知道如何创建、操作和绘制图形。现在我们需要对动力学进行编程。我们的动力学将以一个步骤的方式运作。我们已经概述了每一步会发生什么（如果你不记得了，请参考 "传播动态 "部分）。所以我们知道我们的函数会是什么样子的。
function update_nodes!()
# compute the neighbors of the infected vertices that are susceptible
# update the infected vertices (I -> R)
# infect (S -> I) a subset of susceptible neighbors of infected vertices based on some rule
end
我们可以通过在vertices(G) 循环来计算顶点（节点）的状态。
function get_nodes_state(G)
S = [] # susceptible
I = [] # infected
R = [] # recovered
for v in vertices(G)
state = get_prop(G, v, :state)
if state=="S"
push!(S, v)
elseif state=="I"
push!(I, v)
else
push!(R, v)
end
end
return S, I, R
end
在get_nodes_state 中，我们正在循环浏览G 中的所有顶点，获得它们的道具，并根据它们的:state ，将它们分配到一个列表S 、I 或R 。
首先，update_nodes! 将计算受感染节点的易受感染的邻居。
function update_nodes!(G, I)
# looping through the infected vertices
for v in I
potential_infected_neighbors = []
# looping through their neighbors
for u in neighbors(G, v)
if get_prop(G, u, :state)=="S"
# a susceptible neighbor goes to potential_infected_neighbors
push!(potential_infected_neighbors, u)
end
end
end
end
我们把图G 和之前计算的受感染顶点列表I 传递给update_nodes! 。它将循环浏览I ，计算受感染顶点的所有邻居，如果他们是易受感染的，就把他们加入列表potential_infected_neighbors 。
让我们创建两个函数来更新顶点的状态。
# (I -> R)
function set_node_recovered!(G, node)
set_props!(G, node, Dict(:state=>"R", :color=>"blue"))
end
# (S -> I)
function set_node_infectious!(G, node)
set_props!(G, node, Dict(:state=>"I", :color=>"red"))
end
函数set_nodes_recoreved! 将改变要恢复的顶点的道具。同样，set_node_infectious 将改变顶点的道具为受感染。
如前所述，上一步中被感染的节点在下一步中恢复了。因此，在update_nodes! ，我们已经可以恢复我们已经计算过的被感染的节点。让我们在update_nodes! 中添加set_node_recovered!(G, v) 。这就是它应该有的样子。
function update_nodes!(G, I)
# looping through the infected vertices
for v in I
# updating the infected nodes
set_node_recovered!(G, v)
potential_infected_neighbors = []
# looping through their neighbors
for u in neighbors(G, v)
if get_prop(G, u, :state)=="S"
# a susceptible neighbor goes to potential_infected_neighbors
push!(potential_infected_neighbors, u)
end
end
end
end
我们正在计算受感染的节点，恢复它们并得到它们的易受感染的邻居。
现在我们需要设置一个规则，说明我们应该如何选择下一步将被感染的邻居。让我们保持简单：每个受感染的顶点将感染r 随机邻居。如果顶点的邻居少于r ，它就会感染其所有的邻居。
为了选择r 的随机邻居，让我们使用StatsBase 包中的sample 函数。导入StatsBase 。
using StatsBase
我们可以使用sample(potential_infected_neighbors, r) ，得到一个大小为r 的样本。让我们创建一个列表来存储将要被感染的顶点，并将参数r 传给update_nodes! 。
function update_nodes!(G, I, r)
next_infected_nodes = []
# rest of function here...
end
现在我们可以将potential_infected_neighbors 样本添加到next_infected_nodes 。这就是update_nodes! 的样子。
function update_nodes!(G, I, r)
next_infected_nodes = []
# looping through the infected vertices
for v in I
# updating the infected nodes
set_node_recovered!(G, v)
potential_infected_neighbors = []
# looping through their neighbors
for u in neighbors(G, v)
if get_prop(G, u, :state)=="S"
# a susceptible neighbor goes to potential_infected_neighbors
push!(potential_infected_neighbors, u)
end
end
# getting the samples
l = length(potential_infected_neighbors)
if r <= l
# if there are more than r neighbors, we add the sample to next_infected_nodes
union!(next_infected_nodes, sample(potential_infected_neighbors, r))
else
# if there are less than r neighbors, all the neighbors will be infected
union!(next_infected_nodes, potential_infected_neighbors)
end
end
end
r 如果potential_infected_neighbors 中有超过r 的元素，我们将得到一个大小为potential_infected_neighbors 的样本。如果没有，那么potential_infected_neighbors 中的所有顶点都被感染了。当我们得到这个样本时，我们与next_infected_nodes 进行联合，以防止两个重复的顶点被感染。
然后，我们必须使用next_infected_nodes 来感染。
for v in next_infected_nodes
set_node_infectious!(G, v)
end
在这个片段中，我们通过next_infected_nodes 循环计算，将set_node_infectious! 应用于next_infected_nodes 的每个顶点。现在我们已经有了下一步的顶点和被感染的顶点。
现在，我们有如下的update_nodes! 。
function update_nodes!(G, I, r)
next_infected_nodes = []
# looping through the infected vertices
for v in I
# updating the infected nodes
set_node_recovered!(G, v)
potential_infected_neighbors = []
# looping through their neighbors
for u in neighbors(G, v)
if get_prop(G, u, :state)=="S"
# a susceptible neighbor goes to potential_infected_neighbors
push!(potential_infected_neighbors, u)
end
end
# getting the samples
l = length(potential_infected_neighbors)
if r <= l
# if there are more than r neighbors, we add the sample to next_infected_nodes
union!(next_infected_nodes, sample(potential_infected_neighbors, r))
else
# if there are less than r neighbors, all the neighbors will be infected
union!(next_infected_nodes, potential_infected_neighbors)
end
end
for v in next_infected_nodes
set_node_infectious!(G, v)
end
end
现在让我们创建一个函数来计算下一步并绘制新的图形。
function spreading_plot!(G, r)
S, I, R = get_nodes_state(G)
update_nodes!(G, I, r)
return plot_graph(G)
end
在这个函数中，我们通过运行get_nodes_state 和update_nodes! 来计算下一步的图的状态。然后我们用plot_graph 绘制新的状态。
现在你可以感染一些顶点，为r 选择一个值，然后运行spreading_plot!!试试这个。
set_node_infectious!(G, 1)
r = 1
spreading_plot!(G, r)
我们在G 中感染了顶点1 。我们还设置了r = 1 ，这意味着每个顶点将感染1个邻居。然后我们使用spreading_plot 计算并绘制出图形的下一个状态。
你应该看到像这样的东西。
我们感染了顶点1并运行了传播过程。它有两个邻居，但由于r = 1 ，它只感染了其中一个：图中红色的顶点。顶点1也恢复了，因为在上一个步骤中它被感染了。其余的顶点仍然是易受影响的。
第5步：生成一个动画
如果我们能在一个更大的图形中看到一步步发生的传播过程，那就更好了，对吗？所以，我们现在就来做这个吧
导入这些包。
using Plots, Cairo, Fontconfig, Compose, Printf
创建动画对象。
anim = Animation()
让我们定义一个步骤数t 。这就是我们要运行的次数sreading_plot 。
t = 10
现在我们为我们的动画生成帧--每一帧是一个步骤。
for i in 1:t
# getting plot
p = spreading_plot!(G, r)
output = compose(p, (context(), Compose.rectangle(-10,-10,20,20), fill("white"), Compose.stroke("black")))
# saving the frames
tmpfilename=joinpath(anim.dir, @sprintf("%06d.png",i))
Compose.draw(PNG(tmpfilename),output)
# adding the frames to our animation
push!(anim.frames, tmpfilename)
end
我们在for 循环内运行代码t 次。在这段代码中，我们运行spreading_plot ，在每个迭代中推进模拟的1个时间步。通过compose ，我们在绘图的背景上添加一个白色的矩形。我们为绘图文件创建一个临时路径，并将绘图保存在这个路径上。然后我们在最后一行将其添加到我们的动画anim 的帧中。
为了创建一个gif，我们这样做。
gif(anim, "spreading_infection.gif", fps = 1)
让我们把它变成一个函数!我们只需要把所有的东西都放在函数上make_anim 。
function make_anim(G, r, t)
anim = Animation()
for i in 1:t
p = spreading_plot!(G, r)
output = compose(p, (context(), Compose.rectangle(-10,-10,20,20), fill("white"), Compose.stroke("black")))
tmpfilename=joinpath(anim.dir, @sprintf("%06d.png",i))
Compose.draw(PNG(tmpfilename),output)
push!(anim.frames, tmpfilename)
end
return gif(anim, "spreading_infection.gif", fps = 1)
end
现在让我们创建一个更大的图形，看看我们在本教程中的所有工作是否有效生成一个有100个顶点的Watts-Strogatz图。
# Watts-Strogatz graph with 100 vertices,
# expected degree 4 and edges randomized with probability 0.3
G = MetaGraph(watts_strogatz(100, 4, 0.3))
启动顶点状态，将顶点1设置为感染性。
init_nodes_states!(G)
set_node_infectious!(G, 1)
按照你的意愿设置t 和r 。比如说。
t = 100
r = 2
现在只需运行。
make_anim(G, r, τ)
你应该看到像这样的东西。
看看感染的增长速度有多快!
总结
在本教程中，我们学习了如何使用Graphs.jl和MetaGraphs.jl包在Julia中创建图形并对其进行操作。我们还看到了图是如何成为一个强大的建模工具的，因为它可以被用来模拟对象或代理之间的关系。
利用图的这一特性，我们创建了一个简单的流行病传播模型。我们还学习了如何使用GraphPlot.jl创建图形的可视化。最后，我们把这一切放在一起，创建了一个动画，显示了流行病的传播步骤。
这就是了!这是一个非常简单的模型，但请你花点时间玩一玩，做一些你自己的修改。
这里有一些建议。
- 如果顶点在一个以上的时间步骤中保持感染性会怎样？
- 如果
r更大或更小会怎样？ - 如果接触网时常变化呢？
- 如果被感染的节点没有恢复而是再次变得易受感染呢？
- 如果有些顶点是免疫的（例如接种了疫苗），会怎么样？
