Chinaunix首页 | 论坛 | 博客
  • 博客访问: 3262405
  • 博文数量: 287
  • 博客积分: 0
  • 博客等级: 民兵
  • 技术积分: 2110
  • 用 户 组: 普通用户
  • 注册时间: 2019-10-28 13:40
文章分类

全部博文(287)

文章存档

2022年(60)

2021年(155)

2020年(50)

2019年(22)

我的朋友

分类: 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

# 在下标0endklist列表中根据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(下标为0end 之间)中寻找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()

阅读(295) | 评论(0) | 转发(1) |
给主人留下些什么吧!~~