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()