前言

说实话，已经好久没写博文了，上篇博文的时间都是两个月前。但在这种情况下，粉丝蹭蹭往上涨，hh，这属实是我没想到的，让我着实有种愧疚感。
当然最近没写博文，除了找不到合适的题材外，最重要的是最近有亿点点忙（上课、实验报告、健身、做项目、开会…）剩下的时间都是零零散散，没办法静下心来写东西。

这次主要是趁着要写实验报告，所以顺手将实验报告的思路改写成博文，而作业题目也挺有意思的，是耦合网络信息传播，可以模拟病毒扩散或者信息的扩散，听着是不是有种高大上的感觉？所以我把它写作博文供大家学习参考。（又水了一篇~）

注：本博文代码仅供学习交流，谢绝抄袭！

一、题目介绍

作业题目是一个PPT文件，上来就是一个大标题

第二页就是解题思路（说实话有些突兀，我连题目都不知道是啥就开始叫我一步步做下去，属实有点懵了，这里最好加一个具体的问题，比如用Python编写程序模拟耦合网络的信息传播）

以上便是题目的PPT内容。

二、解题思路

题目内容大致是要我们编写程序模拟耦合网络的信息传播。

解题思路PPT文件其实已经告诉我们了，我们只需按照PPT的步骤一步步进行代码复现即可。

1.构造一个耦合网络

按照PPT的意思，应该是要我们构造一个ER-BA双层网络。

①构造ER网络

生成ER随机网络(算法):
初始: 网络总结点数 N, 没有边；
网络构建过程：每对节点以概率 p 被选择，进行连边，不允许重复连边。

我这里用邻接矩阵（对于Python就是一个二维列表）来存储网络，ER[i][j]==1就说明i和j节点之间有连线，反之无连线。这里用一个全局二维数组变量ER来存储。

根据生成的规则我们编写如下代码来生成ER随机网络（注意这里方法的封装，把这些可变的概率、节点数当做入参，以便后续改变参数），这里我们利用random.random()方法来模拟概率，当随机生成的数字小于概率就生成该边。

# ER随机网络
ER = [[]]
# 生成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

②构造BA网络

BA无标度网络

(1) 初始时网络有 m_0 个节点，可以没有边，也可以像ER网络一样生成一些边；也可以前m_0 个节点是全连通结构。 (2)
每个时间步，一个新节点加入这个网络，并从已存在的网络中选择 m个节点(m≤m_0) 与之相连，选择概率 p_i 与已存在的节点 i的度成正比 p_i=k_i/∑_j▒k_j

(3) 当网络中存在 N 个节点后，停止。

同样BA无标度网络也是以邻接矩阵来存储，这里按照给定的规则进行编码。
考虑到规则比较复杂，可以分为以下几个步骤：

• 初始化BA网络，此时BA网络只有m0个初始节点，边关系随意
• 迭代增加m个节点，每增加一个节点便去根据规则概率寻找m个节点

上述逻辑实现有些复杂，这里我们简化条件来实现题目要求的效果：

• 假设每次增加节点时至少有m个度不为0的节点（因为如果没有m个度不为0的节点，这就意味着我们需要考虑度为0的情况）
• 初始m0个节点组成的网络为全连通网络

由于每次增加节点时选取m个节点操作过于复杂，所以我将其封装为一个方法来增加代码可读性。

以下是代码实现：

# BA无标度网络
BA = [[]]
# 在下标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-BA网络模型

构建出双层网络结构（如ER-ER网卡，ER-BA网络，BA-BA网络），假设层间节点一对一连接，如下图所示。可以理解不同社交平台相互耦合。

这里不需要代码操作，但是我们需要理解，ER节点和BA节点是一一对应的，就像上图展示的那样。
在代码中可以理解为ER网络中下标为i的节点和BA网络中下标为i的节点实际上是连通的。

2.利用SIR模型来模拟信息传播

SIR (Susceptible – Infected - Recoved) model:

SIR传播模型思路:

如果一个 S （健康或者没有接收到信息）状态节点 与一个I （感染或信息传播者）状态节点接触，那么这个S状态的节点被感染的概率为 β .
因此，在t时刻，如果S状态节点有 m 个I状态邻居, 那么下个时间步 t + 1，该节点被感染的概率为 1−(1−β)^m。

t 时刻每个I状态节点在下一时间步，即 t + 1 康复的概率为 γ。

我们可以将上述思路转化为以下步骤：
1.初始一个感染点
2.迭代循环t模拟时间步，每次迭代，每个节点都会根据规则刷新状态

• 如果节点状态为S（正常，未被感染），则它有1−(1−β)^m的概率被感染成I状态
• 如果节点状态为I（被感染），则它有γ的概率康复转变为R状态
• 如果节点状态为R，则状态不会改变（模拟病毒免疫）

3.记录每个时间点的S节点数、I节点数、R节点数

在这里需要注意的是，在考虑周围感染节点时，别忘了考虑ER-BA网络“上下层”的感染情况，因为这里ER网络和BA网络对应的节点时连通的。

这里写代码时也要考虑到入参，同时对于复杂的逻辑进行另外的封装以实现良好的代码可读性，以下是代码实现：

# ER网络感染的情况，S表示为健康，I表示已感染，R表示已康复
ER_SIR = []

BA_SIR = []

# 用于记录每日数据的数组
slist = []
ilist = []
rlist = []
# 对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天，模拟过了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)

3.画图

最后便是根据要求将结果以图表方式可视化展现，这里用了matplotlib这个包，以下是代码实现：

# 画图
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()

三、完整代码

以下是完整代码实现：

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天，模拟过了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()

该代码仅供学习交流，谢绝抄袭！

四、结果分析

输入PPT建议的入参

运行，我们便可以得到一张折线图。

当然，在给出结果之前，我要说明一下初始感染策略——随机感染ER网络中的一个节点，这个前提很重要，会很大程度影响结果。

这运行结果非常有意思,主要为以下两种结果：

1.“感染了大部分节点，最终病毒免疫”

①图表分析

可以看到，在前五天，病毒从最初的一个迅速传染，达到初始节点数的1/2左右，感染数达到顶峰，而康复的人数不断增加，最终在2000节点数的情况下，10天左右有1700个左右的节点感染了病毒并康复，同时还剩300人不到没有感染，此时疫情得到遏止。

多次运行下来，参数不变的情况下，疫情基本都会在感染节点数为1700左右便消失。

②原因猜测

在感染了大多数节点的同时，节点是有几率康复的，而康复后便无法再感染病毒（即R状态无法转换到I状态），所以在感染了大多数节点后，大多数节点陆续康复，在图表中显示的便是R节点数量不断增加。

而为什么还有300个节点不到没有感染呢？
这个大概是因为整个网络的连通度有关，有些节点与其他节点接触不多，在仅有的几次接触中没有被感染，而被感染的节点恰好康复了，所以就存在某些节点没有被感染。

2.“病毒并未传染开便已消失”

①图表分析

这个情况就很有意思了，感染节点数始终为初始节点数1，也就是说病毒就从未传染开。

②原因猜测

在ER随机网络中，由于设置的参数比较小，所以连通性还是比较差的，而恰巧这次初始感染节点是一个和其他节点连通性差的节点，同时它还未来得及传染其他节点时便康复了，此时病毒就已经消失了，在图表中R节点从始至终便为1。

总结

其实我Python挺菜的，写这次作业的时候也是边学边写，写代码的思路很大程度还是照搬Java那套。当然语言是相通的，只是人类操作计算机的工具，所以这个过程不算太难。

这次作业还是蛮有意思的，而运行结果也是出乎我的意料，但仔细分析后也不难推测出原因。
从这次实验中难免联想到现今的疫情，当然现实肯定比题目中的设定复杂，而参数逻辑也会不同，比如病毒会变异，不存在“病毒免疫”的逻辑假设，易感染系数也比设定的要强等等。但还是从一定程度上模拟了流感传播的过程。

注：本文代码仅供学习交流，谢绝抄袭

愿我们以梦为马，不负青春韶华！
与君共勉