用Python制作在地圖上模擬瘟疫擴(kuò)散的Gif圖
受杰森的《Almost Looks Like Work》啟發(fā),我來(lái)展示一些病毒傳播模型。需要注意的是這個(gè)模型并不反映現(xiàn)實(shí)情況,因此不要誤以為是西非可怕的傳染病。相反,它更應(yīng)該被看做是某種虛構(gòu)的僵尸爆發(fā)現(xiàn)象。那么,讓我們進(jìn)入主題。
這就是SIR模型,其中字母S、I和R反映的是在僵尸疫情中,個(gè)體可能處于的不同狀態(tài)。
- S 代表易感群體,即健康個(gè)體中潛在的可能轉(zhuǎn)變的數(shù)量。
- I 代表染病群體,即僵尸數(shù)量。
- R 代表移除量,即因死亡而退出游戲的僵尸數(shù)量,或者感染后又轉(zhuǎn)回人類(lèi)的數(shù)量。但對(duì)與僵尸不存在治愈者,所以我們就不要自我愚弄了(如果要把SIR模型應(yīng)用到流感傳染中,還是有治愈者的)。
- 至于β(beta)和γ(gamma):
- β(beta)表示疾病的傳染性程度,只要被咬就會(huì)感染。
- γ(gamma)表示從僵尸走向死亡的速率,取決于僵尸獵人的平均工作速率,當(dāng)然,這不是一個(gè)完美的模型,請(qǐng)對(duì)我保持耐心。
- S′=?βIS告訴我們健康者變成僵尸的速率,S′是對(duì)時(shí)間的導(dǎo)數(shù)。
- I′=βIS?γI告訴我們感染者是如何增加的,以及行尸進(jìn)入移除態(tài)速率(雙關(guān)語(yǔ))。
- R′=γI只是加上(gamma I),這一項(xiàng)在前面的等式中是負(fù)的。
上面的模型沒(méi)有考慮S/I/R的空間分布,下面來(lái)修正一下!
一種方法是把瑞典和北歐國(guó)家分割成網(wǎng)格,每個(gè)單元可以感染鄰近單元,描述如下:
其中對(duì)于單元,和是它周?chē)乃膫€(gè)單元。(不要因?yàn)閷?duì)角單元而腦疲勞,我們需要我們的大腦不被吃掉)。
初始化一些東東。
import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from matplotlib import rcParams import matplotlib.image as mpimg rcParams['font.family'] = 'serif' rcParams['font.size'] = 16 rcParams['figure.figsize'] = 12, 8 from PIL import Image
適當(dāng)?shù)腷eta和gamma值就能夠摧毀大半江山
beta = 0.010 gamma = 1
還記得導(dǎo)數(shù)的定義么?當(dāng)導(dǎo)數(shù)已知,假設(shè)Δt很小的情況下,經(jīng)過(guò)重新整理,它可以用來(lái)近似預(yù)測(cè)函數(shù)的下一個(gè)取值,我們已經(jīng)聲明過(guò)u′(t)。
初始化一些東東。
import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from matplotlib import rcParams import matplotlib.image as mpimg rcParams['font.family'] = 'serif' rcParams['font.size'] = 16 rcParams['figure.figsize'] = 12, 8 from PIL import Image
適當(dāng)?shù)腷eta和gamma值就能夠摧毀大半江山
beta = 0.010 gamma = 1
還記得導(dǎo)數(shù)的定義么?當(dāng)導(dǎo)數(shù)已知,假設(shè)Δt很小的情況下,經(jīng)過(guò)重新整理,它可以用來(lái)近似預(yù)測(cè)函數(shù)的下一個(gè)取值,我們已經(jīng)聲明過(guò)u′(t)。
這種方法叫做歐拉法,代碼如下:
def euler_step(u, f, dt): return u + dt * f(u)
我們需要函數(shù)f(u)。友好的numpy提供了簡(jiǎn)潔的數(shù)組操作。我可能會(huì)在另一篇文章中回顧它,因?yàn)樗鼈兲珡?qiáng)大了,需要更多的解釋?zhuān)F(xiàn)在這樣就能達(dá)到效果:
def f(u): S = u[0] I = u[1] R = u[2] new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + S[0:-2, 1:-1]*I[0:-2, 1:-1] + S[2:, 1:-1]*I[2:, 1:-1] + S[1:-1, 0:-2]*I[1:-1, 0:-2] + S[1:-1, 2:]*I[1:-1, 2:]), beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + S[0:-2, 1:-1]*I[0:-2, 1:-1] + S[2:, 1:-1]*I[2:, 1:-1] + S[1:-1, 0:-2]*I[1:-1, 0:-2] + S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1], gamma*I[1:-1, 1:-1] ]) padding = np.zeros_like(u) padding[:,1:-1,1:-1] = new padding[0][padding[0] < 0] = 0 padding[0][padding[0] > 255] = 255 padding[1][padding[1] < 0] = 0 padding[1][padding[1] > 255] = 255 padding[2][padding[2] < 0] = 0 padding[2][padding[2] > 255] = 255 return padding
導(dǎo)入北歐國(guó)家的人口密度圖并進(jìn)行下采樣,以便較快地得到結(jié)果
from PIL import Image img = Image.open('popdens2.png') img = img.resize((img.size[0]/2,img.size[1]/2)) img = 255 - np.asarray(img) imgplot = plt.imshow(img) imgplot.set_interpolation('nearest')
北歐國(guó)家的人口密度圖(未包含丹麥)
S矩陣,也就是易感個(gè)體,應(yīng)該近似于人口密度。感染者初始值是0,我們把斯德哥爾摩作為第一感染源。
S_0 = img[:,:,1] I_0 = np.zeros_like(S_0) I_0[309,170] = 1 # patient zero
因?yàn)檫€沒(méi)人死亡,所以把矩陣也置為0.
R_0 = np.zeros_like(S_0)
接著初始化模擬時(shí)長(zhǎng)等。
T = 900 # final time dt = 1 # time increment N = int(T/dt) + 1 # number of time-steps t = np.linspace(0.0, T, N) # time discretization # initialize the array containing the solution for each time-step u = np.empty((N, 3, S_0.shape[0], S_0.shape[1])) u[0][0] = S_0 u[0][1] = I_0 u[0][2] = R_0
我們需要自定義一個(gè)顏色表,這樣才能將感染矩陣顯示在地圖上。
import matplotlib.cm as cm theCM = cm.get_cmap("Reds") theCM._init() alphas = np.abs(np.linspace(0, 1, theCM.N)) theCM._lut[:-3,-1] = alphas
下面坐下來(lái)欣賞吧…
for n in range(N-1): u[n+1] = euler_step(u[n], f, dt)
讓我們?cè)僮鲆幌聢D像渲染,把它做成gif,每個(gè)人都喜歡gifs!
from images2gif import writeGif keyFrames = [] frames = 60.0 for i in range(0, N-1, int(N/frames)): imgplot = plt.imshow(img, vmin=0, vmax=255) imgplot.set_interpolation("nearest") imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM) imgplot.set_interpolation("nearest") filename = "outbreak" + str(i) + ".png" plt.savefig(filename) keyFrames.append(filename) images = [Image.open(fn) for fn in keyFrames] gifFilename = "outbreak.gif" writeGif(gifFilename, images, duration=0.3) plt.clf()
相關(guān)文章
python自動(dòng)結(jié)束mysql慢查詢(xún)會(huì)話的實(shí)例代碼
這篇文章主要介紹了python自動(dòng)結(jié)束mysql慢查詢(xún)會(huì)話,主要涉及到了mysql慢查詢(xún)會(huì)話查詢(xún),定時(shí)任務(wù)的相關(guān)知識(shí),本文通過(guò)實(shí)例代碼給大家介紹的非常詳細(xì),需要的朋友可以參考下2019-10-10在Python中預(yù)先初始化列表內(nèi)容和長(zhǎng)度的實(shí)現(xiàn)
今天小編就為大家分享一篇在Python中預(yù)先初始化列表內(nèi)容和長(zhǎng)度的實(shí)現(xiàn),具有很好的參考價(jià)值,希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2019-11-11Python for Informatics 第11章之正則表達(dá)式(四)
這篇文章主要介紹了Python for Informatics 第11章之正則表達(dá)式(四) 的相關(guān)資料,需要的朋友可以參考下2016-04-04python進(jìn)行數(shù)據(jù)合并concat/merge
這篇文章主要介紹了python進(jìn)行數(shù)據(jù)合并concat/merge,文章圍繞主題展開(kāi)詳細(xì)的內(nèi)容介紹,具有一定的參考價(jià)值,感興趣的小伙伴可以參考一下2022-09-09PIL.Image.open和cv2.imread的比較與相互轉(zhuǎn)換的方法
這篇文章主要介紹了PIL.Image.open和cv2.imread的比較與相互轉(zhuǎn)換的方法,文中通過(guò)示例代碼介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2020-06-06