分类： Python/Ruby
2021-12-16 12:08:10
import random
import matplotlib.pyplot as plt
# ER随机网络
ER = [[]],
# BA无标度网络
BA = [[]]
# ER网络感染的情况，S表示为健康，I表示已感染，R表示已康复
ER_SIR = []
BA_SIR = []
# 用于记录每日数据的数组
slist = []
ilist = []
rlist = []
# 生成ER随机网络,n表示点个数，p表示生成边的概率
def generateER(n, p):
global ER
# 避免浅拷贝
ER = [([0] * n) for i in range(n)]
for i in range(0, n):
for j in range(i + 1, n):
# 随机生成的数
x = random.random()
# 概率生成边
if x < p:
ER[i][j] = 1
ER[j][i] = 1
# 在下标0到end的klist列表中根据k的权重来寻找一个下标返回,不包括exclude_list中的节点
def findOne(klist, end, exclude_list):
flag = True
k_sum = 0
for i in range(0, end):
k_sum += klist[i]
result = end - 1
while flag:
x = random.random() * k_sum
pre = 0
for i in range(0, end):
pre += klist[i]
if pre > x:
result = i
break
flag = False
# 检查结果是否在exclude_list中，如果是则重来，否则返回结果
for i in exclude_list:
if i == result:
flag = True
break
return result
# 根据BA无标度网络的规则在klist(下标为0到end 之间)中寻找m个数
# 为了简化逻辑，这里参与竞争（k！=0）的节点数必定大于m
def find(k_list, m, end):
# 结果数组
result = []
# 初始化
for i in range(0, m):
j = findOne(k_list, end, result)
# 增加结果
result.append(j)
return result
# 生成BA无标度网络,n表示点个数，m0表示初始节点数
def generateBA(n, m0, m):
global BA
# 初始BA无标度网络(利用生成器创建避免浅拷贝)
BA = [([0] * n) for i in range(n)]
# 初始化前m0个节点都为连通网络
for i in range(0, m0):
for j in range(i + 1, m0):
BA[i][j] = 1
BA[j][i] = 1
# 初始度数组
k_list = [m0 - 1] * m0 + [0] * (n - m)
# 遍历后面的节点，模拟加入
for i in range(m0, n):
# 选出m个节点
nodes = find(k_list, m, i)
for j in nodes:
BA[i][j] = 1
BA[j][i] = 1
# 更新度数组
k_list[i] += 1
k_list[j] += 1
# 对ER网络的i节点进行处理，i表示节点下标，t表示天数，b表示感染系数，y表示康复系数
def er_sir_one(i, t, b, y):
global ER, BA, ER_SIR, BA_SIR, slist, ilist, rlist
if ER_SIR[i] == 'S':
# 开始统计周围节点感染的数量
inum = 0
for x in range(len(ER_SIR)):
if (not x == i) and ER[i][x] == 1 and ER_SIR[x] == 'I':
inum += 1
# 因为是双层网络，所以BA也要考虑
if BA_SIR[i] == 'I':
inum += 1
ilist[t] += 1
p = random.random()
# 有1-(1-b)^inum概率感染
if p < 1 - (1 - b) ** inum:
ER_SIR[i] = 'I'
ilist[t] += 1
return
slist[t] += 1
elif ER_SIR[i] == 'I':
p = random.random()
# 有y的几率康复
if p < y:
ER_SIR[i] = 'R'
rlist[t] += 1
return
ilist[t] += 1
else:
rlist[t] += 1
# 对BA网络的i节点进行处理，i表示节点下标，t表示天数，b表示感染系数，y表示康复系数
def ba_sir_one(i, t, b, y):
global ER, BA, ER_SIR, BA_SIR, slist, ilist, rlist
if BA_SIR[i] == 'S':
# 开始统计周围节点感染的数量
inum = 0
for x in range(len(BA_SIR)):
if (not x == i) 外汇跟单gendan5.comand BA[i][x] == 1 and BA_SIR[x] == 'I':
inum += 1
# 因为是双层网络，所以ER也要考虑
if ER_SIR[i] == 'I':
inum += 1
p = random.random()
# 有1-(1-b)^inum概率感染
if p < 1 - (1 - b) ** inum:
BA_SIR[i] = 'I'
ilist[t] += 1
return
slist[t] += 1
# 有y的几率康复
elif BA_SIR[i] == 'I':
p = random.random()
if p < y:
BA_SIR[i] = 'R'
rlist[t] += 1
return
ilist[t] += 1
else:
rlist[t]+=1
def sir(b, y, t):
global ER_SIR, BA_SIR, slist, ilist, rlist
n = len(ER[0])
# 初始化ER_SIR,BA_SIR
ER_SIR = ['S'] * n
BA_SIR = ['S'] * n
# 初始化每日统计数据的数组
slist = [0] * t
ilist = [0] * t
rlist = [0] * t
# 随机感染ER网络中的一个节点
x = random.randint(0, n - 1)
ER_SIR[x] = 'I'
# 遍历t天，模拟过了t天
for day in range(t):
# 遍历所有节点
for node in range(n):
# 对双层网络状态进行一遍刷新
er_sir_one(node, day, b, y)
ba_sir_one(node, day, b, y)
# 画图
plt.plot(list(range(t)), slist, linewidth=4,label=u'S')
plt.plot(list(range(t)), ilist, linewidth=4,label=u'I')
plt.plot(list(range(t)), rlist, linewidth=4,label=u'R')
plt.legend() # 让图例生效
plt.xlabel(u"t") # X轴标签
plt.ylabel("N(t)") # Y轴标签
plt.title("SIR model result") # 标题
plt.show()
def main():
generateER(1000, 0.006)
generateBA(1000, 3, 3)
sir(0.2, 0.5, 50)
if __name__ == '__main__':
main()
