前言

說實話，已經好久沒寫博文了，上篇博文的時間都是兩個月前，但在這種情況下，粉絲蹭蹭往上漲，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那套，當然語言是相通的，只是人類操作計算機的工具，所以這個程序不算太難，

這次作業還是蠻有意思的，而運行結果也是出乎我的意料，但仔細分析后也不難推測出原因，
從這次實驗中難免聯想到現今的疫情，當然現實肯定比題目中的設定復雜，而引數邏輯也會不同，比如病毒會變異，不存在“病毒免疫”的邏輯假設，易感染系數也比設定的要強等等，但還是從一定程度上模擬了流感傳播的程序，

注：本文代碼僅供學習交流，謝絕抄襲

愿我們以夢為馬，不負青春韶華！
與君共勉