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)