SI模型
import numpy as np import matplotlib.pyplot as plt import smallworld as sw #邻接矩阵 a = sw.a # 感染率 beta = sw.beta #初始患者 origin = sw.origin def si_(a, beta, origin): #总人数 n = a.shape[0] #控制符 judge = 1 #未感染人群 s = np.arange(n) s = np.delete(s, origin) #上一轮或原先感染人群 i = origin #感染者记录 r = [] #感染时间记录 t = [] #感染人数记录 speed = [] #计时器 h = 0 #总感染者 infected = i #开始传播 while judge == 1: temp_i = [] #感染人数 m = len(infected) for j in range(0, m): node = int(infected[j]) asd = [] for k in range(0, n): if a[node, k] == 1: asd.append(k) #感染者邻接的未感染者 asd2 = np.intersect1d(asd, s) #随机生成被感染率 num = np.random.rand(len(asd2)) - beta #该感染者传播的新感染者 asd_final =[] for k in range(0, len(asd2)): if num[k] <= 0: asd_final.append(asd2[k]) #这一轮总的新感染者 temp_i = np.union1d(temp_i, asd_final) s = np.setdiff1d(s, asd_final) if len(i) > 0: for k in range(0, len(i)): r.append(i[k]) t.append(h) #新一轮感染人群 i = temp_i infected = np.union1d(infected, i) #所有人都被感染则跳出循环 if len(s) == 0: judge = 0 if len(i) > 0: for k in range(0, len(i)): r.append(i[k]) t.append(h) speed.append(len(r)) h = h+1 coverage = r for j in range(h, 2*h): speed.append(speed[j-1]) print(speed) plt.plot(speed, "-ro", label='Infectious') plt.title('SI') plt.legend() plt.show() si_(a, beta, origin)
效果图：
SIR模型
import numpy as np import matplotlib.pyplot as plt import smallworld as sw #邻接矩阵 a = sw.a #感染率 beta = sw.beta #复原率 gama = sw.gama #初始患者 origin = sw.origin def sir_(a, origin, beta, gama): #总人数 n = a.shape[0] #控制符 judge = 1 #未感染人群 s = np.arange(n) s = np.delete(s, origin) # 未感染人数记录 snum = list() snum.append(len(s)) #感染人群 i = origin #感染人数记录 inum = list() inum.append(len(i)) #复原人群 r = [] #复原人数记录 rnum = list() rnum.append(0) #计时器 h = 1 #开始传播 while judge ==1: #新一轮感染者 temp_i = [] m = len(i) #感染 for j in range(0, m): #感染者 node = int(i[j]) asd = [] for k in range(0, n): if a[node, k] == 1: asd.append(k) #未被感染者且不为复原者 asd2 = np.intersect1d(asd, s) asd2 = np.setdiff1d(asd2, r) num = np.random.rand(len(asd2)) - beta #新感染者 asd_final = [] for k in range(0, len(asd2)): if num[k] <= 0: asd_final.append(asd2[k]) temp_i = np.union1d(temp_i, asd_final) s = np.setdiff1d(s, asd_final) #复原 num = np.random.rand(m) - gama asd = [] asd_final = [] for k in range(0, m): if num[k] <= 0: asd.append(k) asd_final.append(i[k]) #原感染者部分复原 r = np.union1d(r, asd_final) i = np.delete(i, asd) #这一轮后总感染者 i = np.union1d(i, temp_i) snum.append(len(s)) inum.append(len(i)) rnum.append(len(r)) h = h + 1 if len(i) == 0: judge = 0 for j in range(h, h*2): snum.append(snum[j-1]) inum.append(inum[j - 1]) rnum.append(rnum[j - 1]) plt.plot(snum, "-bo", label='Susceptibles') plt.plot(inum, "-ro", label='Infectious') plt.plot(rnum, "-go", label='Recovereds') plt.title('SIR') plt.legend() plt.show() sir_(a, origin, beta, gama)
效果图：
SIS模型
import numpy as np import matplotlib.pyplot as plt import BA import random #邻接矩阵 a = BA.a #感染率 beta = 0.6 #复原率 gama = 0.3 #初始感染者 origin = random.sample(range(0, a.shape[0]), 5) def sis_(a, origin, beta, gama): #总人数 n = a.shape[0] #未感染人群 s = np.arange(n) s = np.delete(s, origin) #感染人群 i = origin #感染人数记录 speed = [] speed.append(len(i)) #计时器 h = 1 while h < 70: temp_i = [] m = len(i) for j in range(0, m): node = int(i[j]) asd = [] for k in range(0, n): if a[node, k] == 1: asd.append(k) asd2 = np.intersect1d(asd, s) num = np.random.rand(len(asd2)) - beta asd_final = [] for k in range(0, len(asd2)): if num[k] <= 0: asd_final.append(asd2[k]) temp_i = np.union1d(temp_i, asd_final) s = np.setdiff1d(s, asd_final) num = np.random.rand(m) - gama asd = [] asd_final = [] for j in range(0, m): if num[j] <= 0: asd.append(j) asd_final.append(i[j]) s = np.union1d(s, asd_final) i = np.delete(i, asd) i = np.union1d(i, temp_i) speed.append(len(i)) h = h + 1 snum = [] for j in range(0, h): snum.append(n - speed[j]) plt.plot(speed, "-ro", label='Infectious') plt.plot(snum, "-bo", label='Susceptibles') plt.title('SIS') plt.legend() plt.show() sis_(a, origin, beta, gama)
效果图：
LT模型
import numpy as np import smallworld as sw import networkx as nx import matplotlib.pyplot as plt #邻接矩阵 a = sw.a #节点度数, 1/b是其他节点对该节点的影响力 b = sw.b #节点阀值 beta = sw.beta #原激活节点 origin = sw.origin #超过beta（如50%）的邻接节点处于激活状态，该节点才会进入激活状态 def lt_(a, b, origin, beta): #节点数 n = a.shape[0] #控制符 judge = 1 #未激活节点 s = np.arange(n) s = np.delete(s, origin) #激活节点 i = origin while judge == 1: #该轮激活节点 temp_i = [] #激活节点个数 m = len(i) for j in range(0, m): node = int(i[j]) asd = [] for k in range(0, n): if a[node, k] == 1: asd.append(k) #找到相邻的未激活节点 asd2 = np.intersect1d(asd, s) asd_final = [] for k in range(0,len(asd2)): num = 0 #该未激活节点相邻的激活节点个数 for t in range(0, m): if a[int(i[t]), asd2[k]] == 1: num = num + 1 if 1 / b[asd2[k]] * num >= beta: asd_final.append(asd2[k]) temp_i = np.union1d(temp_i, asd_final) s = np.setdiff1d(s, asd_final) #将新激活节点合并到原激活节点中 i = np.union1d(i, temp_i) #如果该轮没有新激活节点，那之后都不会再有，跳出循环 if len(temp_i) == 0: judge = 0 #输出新的网络状况 color = [] for j in range(0, n): color.append('b') for j in range(0, len(i)): color[int(i[j])] = 'r' g = nx.from_numpy_matrix(a) nx.draw(g, with_labels=True, node_color=color) plt.show() lt_(a, b, origin, beta)
效果图：
数据生成代码（WS）
import numpy as np import random import networkx as nx import matplotlib.pyplot as plt n = 50 k = 5 p = 0.1 #建立小世界网 def sw_(n, k, p): m = np.zeros((n, n), dtype=int) #节点度数 d = np.zeros(n, dtype=int) #构建环形网络 for i in range(0, n): for j in range(0, n): if j > i and j <= i+k: m[i, j] = 1 m[j, i] = 1 d[i] = d[i] + 1 d[j] = d[j] + 1 if i+k >= n and j <= i+k-n: m[i, j] = 1 m[j, i] = 1 d[i] = d[i] + 1 d[j] = d[j] + 1 #讲规则网络转换为随机网络 for i in range(0, n): for j in range(0, n): if m[i, j] == 1: rand_num1 = np.random.rand(1) if rand_num1 <= p: rand_num2 = np.random.rand(1) temp = np.random.permutation(np.arange(n)) if rand_num2 < 0.5: for l in range(0, n): if i != temp[l] and m[i, temp[l]] == 0: m[i, temp[l]] = 1 m[temp[l], i] = 1 d[j] = d[j] - 1 d[temp[l]] = d[temp[l]] + 1 break else: for l in range(0, n): if j != temp[l] and m[j, temp[l]] == 0: m[j, temp[l]] = 1 m[temp[l], j] = 1 d[i] = d[i] - 1 d[temp[l]] = d[temp[l]] + 1 break m[i, j] = 0 m[j, i] = 0 return m, d a, b = sw_(n, k, p) #LT预设 #预设阀值 beta = 0.5 #初始激活节点 q = np.random.randint(10, 15) origin = random.sample(range(0, n), q)