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) and 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 天，[利率期货](https://www.gendan5.com/ff/if.html) 模仿过了 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()