您好,登錄后才能下訂單哦!
病毒擴(kuò)散仿真程序,用 python 也可以。
概述
事情是這樣的,B 站 UP 主 @ele 實(shí)驗(yàn)室,寫了一個(gè)簡單的疫情傳播仿真程序,告訴大家在家待著的重要性,視頻相信大家都看過了,并且 UP 主也放出了源碼。
因?yàn)槭?Java 開發(fā)的,所以開始我并沒有多加關(guān)注。后來看到有人解析代碼,發(fā)現(xiàn)我也能看懂,然后就琢磨用 Python 應(yīng)該怎么實(shí)現(xiàn)。
Java 版程序淺析
一個(gè)人就是 1 個(gè)(x, y)坐標(biāo)點(diǎn),并且每個(gè)人有一個(gè)狀態(tài)。
public class Person extends Point { private int state = State.NORMAL; }
在每一輪的迭代中,遍歷每個(gè)人,每個(gè)人根據(jù)自身的狀態(tài),做出一定的動(dòng)作,包括:
這些動(dòng)作的具體變更,取決于定義的各種系數(shù)。
一輪迭代完成,打印這些點(diǎn),不同的狀態(tài)對應(yīng)不同的顏色。
繪圖部分直接使用的 Java 繪圖類 Graphics。
Python 版思路
如果我們想用 Python 實(shí)現(xiàn)應(yīng)該怎么做呢?
如果完全復(fù)刻 Java 版本,則每次迭代需遍歷所有人,并計(jì)算和他人距離,這就是 N^2 次計(jì)算。如果是 1000 個(gè)人,就需要循環(huán) 1 百萬次。這個(gè) Python 的性能肯定捉急。
不過 Python 有 numpy ,可以快速的操作數(shù)組。結(jié)合 matplotlib 則可以畫出圖形。
import numpy as np import matplotlib.pyplot as plt
如何模擬人群
為了減少函數(shù)之間互相傳參和使用全局變量,我們也來定義一個(gè)類:
class People(object): def __init__(self, count=1000, first_infected_count=3): self.count = count self.first_infected_count = first_infected_count self.init()
所有人的坐標(biāo)數(shù)據(jù)就是 N 行 2 列的數(shù)組,同時(shí)伴隨一定的狀態(tài):
def init(self): self._people = np.random.normal(0, 100, (self.count, 2)) self.reset()
狀態(tài)值和計(jì)時(shí)器也都是數(shù)組,同時(shí)每次隨機(jī)選取指定數(shù)量的人感染:
def reset(self): self._round = 0 self._status = np.array([0] * self.count) self._timer = np.array([0] * self.count) self.random_people_state(self.first_infected_count, 1)
這里關(guān)鍵的一點(diǎn)是,輔助數(shù)組的大小和人數(shù)保持一致,這樣就能形成一一對應(yīng)的關(guān)系。
狀態(tài)發(fā)生變化的人才順帶記錄時(shí)間:
def random_people_state(self, num, state=1): """隨機(jī)挑選人設(shè)置狀態(tài) """ assert self.count > num # TODO:極端情況下會(huì)出現(xiàn)無限循環(huán) n = 0 while n < num: i = np.random.randint(0, self.count) if self._status[i] == state: continue else: self.set_state(i, state) n += 1 def set_state(self, i, state): self._status[i] = state # 記錄狀態(tài)改變的時(shí)間 self._timer[i] = self._round
通過狀態(tài)值,就可以過濾出人群,每個(gè)人群都是 people 的切片視圖。這里 numpy 的功能相當(dāng)強(qiáng)大,只需要非常簡潔的語法即可實(shí)現(xiàn):
@property def healthy(self): return self._people[self._status == 0] @property def infected(self): return self._people[self._status == 1]
按照既定的思路,我們先來定義每輪迭代要做的動(dòng)作:
def update(self): """每一次迭代更新""" self.change_state() self.affect() self.move() self._round += 1 self.report()
順序和開始分析的略有差異,其實(shí)并不是十分重要,調(diào)換它們的順序也是可以的。
如何改變狀態(tài)
這一步就是更新狀態(tài)數(shù)組 self._status 和 計(jì)時(shí)器數(shù)組 self._timer:
def change_state(self): dt = self._round - self._timer # 必須先更新時(shí)鐘再更新狀態(tài) d = np.random.randint(3, 5) self._timer[(self._status == 1) & ((dt == d) | (dt > 14))] = self._round self._status[(self._status == 1) & ((dt == d) | (dt > 14))] += 1
仍然是通過切片過濾出要更改的目標(biāo),然后全部更新。
這里具體的實(shí)現(xiàn)我寫的非常簡單,沒有引入太多的變量:
在一定周期內(nèi)的 感染者(infected),狀態(tài)置為 確診(confirmed)。 我這里簡單假設(shè)了確診者就被醫(yī)院收治,所以失去了繼續(xù)感染他人的機(jī)會(huì)(見下面)。如果要搞復(fù)雜點(diǎn),可以引入病床,治愈,死亡等狀態(tài)。
如何影響他人
影響別人是整個(gè)程序的性能瓶頸,因?yàn)樾枰?jì)算每個(gè)人之間的距離。
這里繼續(xù)做了簡化,只處理感染者:
def infect_possible(self, x=0., safe_distance=3.0): """按概率感染接近的健康人 x 的取值參考正態(tài)分布概率表,x=0 時(shí)感染概率是 50% """ for inf in self.infected: dm = (self._people - inf) ** 2 d = dm.sum(axis=1) ** 0.5 sorted_index = d.argsort() for i in sorted_index: if d[i] >= safe_distance: break # 超出范圍,不用管了 if self._status[i] > 0: continue if np.random.normal() > x: continue self._status[i] = 1 # 記錄狀態(tài)改變的時(shí)間 self._timer[i] = self._round
可以看到,距離的計(jì)算仍然是通過 numpy 的矩陣操作。但是需要對每一個(gè)感染者單獨(dú)計(jì)算,所以如果感染者較多,python 的處理效率感人。
如何移動(dòng)
_people 是一個(gè)坐標(biāo)矩陣,只要生成移動(dòng)距離矩陣 dt,然后它相加即可。我們可以設(shè)置一個(gè)可移動(dòng)的范圍 width,把移動(dòng)距離控制在一定范圍內(nèi)。
def move(self, width=1, x=.0): movement = self.random_movement(width=width) # 限定特定狀態(tài)的人員移動(dòng) switch = self.random_switch(x=x) movement[switch == 0] = 0 self._people = self._people + movement
這里還需要增加一個(gè)控制移動(dòng)意向的選項(xiàng),仍然是利用了正態(tài)分布概率??紤]到這種場景有可能會(huì)重用,所以特地把這個(gè)方法提取了出來,生成一個(gè)只包含 0 1 的數(shù)組充當(dāng)開關(guān)。
def random_switch(self, x=0.): """隨機(jī)生成開關(guān),0 - 關(guān),1 - 開 x 大致取值范圍 -1.99 - 1.99; 對應(yīng)正態(tài)分布的概率, 取值 0 的時(shí)候?qū)?yīng)概率是 50% :param x: 控制開關(guān)比例 :return: """ normal = np.random.normal(0, 1, self.count) switch = np.where(normal < x, 1, 0) return switch
輸出結(jié)果
有了一切數(shù)據(jù)和變化之后,接下來最重要的事情自然就是圖形化顯示結(jié)果了。直接使用 matplotlib 的散點(diǎn)圖就可以了:
def report(self): plt.cla() # plt.grid(False) p1 = plt.scatter(self.healthy[:, 0], self.healthy[:, 1], s=1) p2 = plt.scatter(self.infected[:, 0], self.infected[:, 1], s=1, c='pink') p3 = plt.scatter(self.confirmed[:, 0], self.confirmed[:, 1], s=1, c='red') plt.legend([p1, p2, p3], ['healthy', 'infected', 'confirmed'], loc='upper right', scatterpoints=1) t = "Round: %s, Healthy: %s, Infected: %s, Confirmed: %s" % \ (self._round, len(self.healthy), len(self.infected), len(self.confirmed)) plt.text(-200, 400, t, ha='left', wrap=True)
實(shí)際效果
啟動(dòng)。
if __name__ == '__main__': np.random.seed(0) plt.figure(figsize=(16, 16), dpi=100) plt.ion() p = People(5000, 3) for i in range(100): p.update() p.report() plt.pause(.1) plt.pause(3)
因?yàn)檫@個(gè)小 demo 主要是個(gè)人用來練手,目前一些參數(shù)沒有完全抽出來。有需要的只能直接改源碼。
后記
從多次實(shí)驗(yàn)的結(jié)果,通過調(diào)整人員的流動(dòng)意愿,流動(dòng)距離等因素,是可以得到直觀的結(jié)論的。
本人也是初次使用 numpy 和 matplotlib,現(xiàn)學(xué)現(xiàn)賣,若有使用不當(dāng)之處請指正。其中的概率參數(shù)設(shè)置 基本沒有科學(xué)依據(jù),僅供 Python 愛好者參考。
總得來說,用 numpy 來模擬病毒感染情況應(yīng)該是能行得通的。但是其中的影響因子還需要仔細(xì)設(shè)計(jì)。性能也是需要考量的問題。
源碼地址
總結(jié)
以上所述是小編給大家介紹的Python 寫了個(gè)新型冠狀病毒疫情傳播模擬程序,希望對大家有所幫助,也非常感謝大家對億速云網(wǎng)站的支持!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場,如果涉及侵權(quán)請聯(lián)系站長郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。