问题描述
我想在 NetworkX 图上模拟和可视化 Sirs 模型(其中恢复/删除的节点可能再次变得易受影响)。每个节点也作为代理运行,并且可以在每个时间步选择以概率 p 进行隔离,从而在该时间步内不会被感染。
我已经用 EoN.Gillespie_simple_contagion 模拟了一个 Sirs 模型,我试图了解我是否可以修改这些方法以包含代理行为,或者我是否需要编写一个定制的方法。
我已经看到可以修改一些 EoN 方法:Modified SIR model
这是我使用的代码:
import networkx as nx
import matplotlib.pyplot as plt
import EoN
from collections import defaultdict
# parameters required for the Sirs model
a = 0.1
b = 0.01
y = 0.001
d = 0.001
# Simple contagions
# the below is based on an example of a SEIR disease (there is an exposed state before becoming infectIoUs)
# from https://arxiv.org/pdf/2001.02436.pdf
Gnp = nx.gnp_random_graph(500,0.005)
H = nx.DiGraph() #For the spontaneous transitions
H.add_edge('I','R',rate = b) # an infected node can be recovered/removed
H.add_edge('I','S',rate = y) # an infected node can become susceptible again
H.add_edge('R',rate = d) # a recovered node can become suscepticle again
J = nx.DiGraph() #for the induced transitions
J.add_edge(('I','S'),('I','I'),rate = a) # a susceptible node can become infected from a neighbour
IC = defaultdict(lambda: 'S')
# set all statuses except one to susceptible. only one node shall begin infected
for node in range(500):
IC[node] = 'S'
IC[0] = 'I'
return_statuses = ('S','I','R')
print('doing Gillespie simulation')
t,S,I,R = EoN.Gillespie_simple_contagion(Gnp,H,J,IC,return_statuses,tmax = 500)
print('done with simulation,Now plotting')
plt.plot(t,label = 'Susceptible')
plt.plot(t,label = 'Infected')
plt.plot(t,R,label = 'Recovered')
plt.xlabel('$t$')
plt.ylabel('Number of nodes')
plt.legend()
plt.show()
解决方法
为此,我们将引入一个新的节点状态 'X'
来表示易受影响的隔离节点。我在 'S'
中添加了从 'X'
到 'X'
的边以及从 'S'
到 'H'
的另一个边。我已将 'X'
添加到 return_statuses
,并绘制了它。
为了简化您的代码,我删除了最初将所有内容分配为 'S'
的循环。因为它使用了 default_dict
,所以它都隐含地以 'S'
开头。
您提到下一个“时间步长”并以概率 p
移动。这不是这段代码中会发生的事情。该代码在连续时间内工作。所以它需要一个速率。它没有具体说明一个节点保持隔离的确切时间,而是给出了它随时离开的可能性的概念。如果这对您没有意义,您应该阅读一下指数分布。
虽然我没有在这里使用它,但您可能想知道 EoN 附带的工具允许您animate the simulation。
import networkx as nx
import matplotlib.pyplot as plt
import EoN
from collections import defaultdict
a = 0.1
b = 0.01
y = 0.001
d = 0.001
to_isolation_rate = 0.05
from_isolation_rate = 1
# Simple contagions
# the below is based on an example of a SEIR disease (there is an exposed state before becoming infectious)
# from https://arxiv.org/pdf/2001.02436.pdf
Gnp = nx.gnp_random_graph(500,0.005)
H = nx.DiGraph() #For the spontaneous transitions
H.add_edge('I','R',rate = b) # an infected node can be recovered/removed
H.add_edge('I','S',rate = y) # an infected node can become susceptible again
H.add_edge('R',rate = d) # a recovered node can become suscepticle again
H.add_edge('S','X',rate = to_isolation_rate)
H.add_edge('X',rate = from_isolation_rate)
J = nx.DiGraph() #for the induced transitions
J.add_edge(('I','S'),('I','I'),rate = a) # a susceptible node can become infected from a neighbour
IC = defaultdict(lambda: 'S')
IC[0] = 'I'
return_statuses = ('S','I','X')
print('doing Gillespie simulation')
t,S,I,R,X = EoN.Gillespie_simple_contagion(Gnp,H,J,IC,return_statuses,tmax = 500)
print('done with simulation,now plotting')
plt.plot(t,label = 'Susceptible')
plt.plot(t,label = 'Infected')
plt.plot(t,label = 'Recovered')
plt.plot(t,X,label = 'Isolating')
plt.xlabel('$t$')
plt.ylabel('Number of nodes')
plt.legend()
plt.show()
,
可以为所需的功能编写定制的方法:
def SIRS_simulation_with_quarantine(graph,params,tmax):
a,b,y,d,q = params
t=0
graph_current_timestep = graph.copy()
graph_previous_timestep = graph.copy()
Infected_nodes = [x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='I']
S = []
I = []
R = []
T = []
while(t<tmax):
for node in graph.nodes():
# initialise random probablilty for that node
p = np.random.rand()
# determine neighbours for that node
neighbors = [n for n in graph_previous_timestep.neighbors(node)]
# depending on the state,check for any state changes and update node attributes
if graph_previous_timestep.nodes[node]['state'][0] == 'S':
if (q > np.random.rand()): # if the node didnt quarantine...
if len(set.intersection(set(neighbors),set(Infected_nodes)))>0: # if any neighbors are infected...
if (p < a):
attrs = {node: {'state': 'I'}}
nx.set_node_attributes(graph_current_timestep,attrs)
elif graph_previous_timestep.nodes[node]['state'][0] == 'I':
if (p < b):
attrs = {node: {'state': 'R'}}
nx.set_node_attributes(graph_current_timestep,attrs)
elif (p < y):
attrs = {node: {'state': 'S'}}
nx.set_node_attributes(graph_current_timestep,attrs)
elif graph_previous_timestep.nodes[node]['state'][0] == 'R':
if (p < d):
attrs = {node: {'state': 'S'}}
nx.set_node_attributes(graph_current_timestep,attrs)
# after the changes to states have occured for all nodes,record number of nodes at each state,increment timestep
Infected_nodes = [x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='I']
S.append(len([x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='S']))
I.append(len(Infected_nodes))
R.append(len([x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='R']))
T.append(t)
graph_previous_timestep = graph_current_timestep.copy()
t = t+1
return T,R