用Python制作在地圖上模擬瘟疫擴散的Gif圖
受杰森的《Almost Looks Like Work》啟發(fā)，我來展示一些病毒傳播模型。需要注意的是這個模型并不反映現(xiàn)實情況，因此不要誤以為是西非可怕的傳染病。相反，它更應該被看做是某種虛構(gòu)的僵尸爆發(fā)現(xiàn)象。那么，讓我們進入主題。
這就是SIR模型，其中字母S、I和R反映的是在僵尸疫情中，個體可能處于的不同狀態(tài)。
- S 代表易感群體，即健康個體中潛在的可能轉(zhuǎn)變的數(shù)量。
- I 代表染病群體，即僵尸數(shù)量。
- R 代表移除量，即因死亡而退出游戲的僵尸數(shù)量，或者感染后又轉(zhuǎn)回人類的數(shù)量。但對與僵尸不存在治愈者，所以我們就不要自我愚弄了（如果要把SIR模型應用到流感傳染中，還是有治愈者的）。
- 至于β(beta)和γ(gamma):
- β(beta)表示疾病的傳染性程度，只要被咬就會感染。
- γ(gamma)表示從僵尸走向死亡的速率，取決于僵尸獵人的平均工作速率，當然，這不是一個完美的模型，請對我保持耐心。
- S′=?βIS告訴我們健康者變成僵尸的速率，S′是對時間的導數(shù)。
- I′=βIS?γI告訴我們感染者是如何增加的，以及行尸進入移除態(tài)速率（雙關語）。
- R′=γI只是加上(gamma I)，這一項在前面的等式中是負的。
上面的模型沒有考慮S/I/R的空間分布，下面來修正一下！
一種方法是把瑞典和北歐國家分割成網(wǎng)格，每個單元可以感染鄰近單元，描述如下：
其中對于單元，和是它周圍的四個單元。（不要因為對角單元而腦疲勞，我們需要我們的大腦不被吃掉）。
初始化一些東東。
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
適當?shù)腷eta和gamma值就能夠摧毀大半江山
beta = 0.010 gamma = 1
還記得導數(shù)的定義么？當導數(shù)已知，假設Δt很小的情況下，經(jīng)過重新整理，它可以用來近似預測函數(shù)的下一個取值，我們已經(jīng)聲明過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
適當?shù)腷eta和gamma值就能夠摧毀大半江山
beta = 0.010 gamma = 1
還記得導數(shù)的定義么？當導數(shù)已知，假設Δt很小的情況下，經(jīng)過重新整理，它可以用來近似預測函數(shù)的下一個取值，我們已經(jīng)聲明過u′(t)。
這種方法叫做歐拉法，代碼如下：
def euler_step(u, f, dt): return u + dt * f(u)
我們需要函數(shù)f(u)。友好的numpy提供了簡潔的數(shù)組操作。我可能會在另一篇文章中回顧它，因為它們太強大了，需要更多的解釋，但現(xiàn)在這樣就能達到效果：
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
導入北歐國家的人口密度圖并進行下采樣，以便較快地得到結(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')
北歐國家的人口密度圖（未包含丹麥）
S矩陣，也就是易感個體，應該近似于人口密度。感染者初始值是0，我們把斯德哥爾摩作為第一感染源。
S_0 = img[:,:,1] I_0 = np.zeros_like(S_0) I_0[309,170] = 1 # patient zero
因為還沒人死亡，所以把矩陣也置為0.
R_0 = np.zeros_like(S_0)
接著初始化模擬時長等。
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
我們需要自定義一個顏色表，這樣才能將感染矩陣顯示在地圖上。
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
下面坐下來欣賞吧…
for n in range(N-1): u[n+1] = euler_step(u[n], f, dt)
讓我們再做一下圖像渲染，把它做成gif，每個人都喜歡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()
相關文章
今天小編就為大家分享一篇python繪制直方圖和密度圖的實例，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2019-07-07
python自動結(jié)束mysql慢查詢會話的實例代碼
這篇文章主要介紹了python自動結(jié)束mysql慢查詢會話,主要涉及到了mysql慢查詢會話查詢，定時任務的相關知識，本文通過實例代碼給大家介紹的非常詳細，需要的朋友可以參考下2019-10-10
在Python中預先初始化列表內(nèi)容和長度的實現(xiàn)
今天小編就為大家分享一篇在Python中預先初始化列表內(nèi)容和長度的實現(xiàn)，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2019-11-11
這篇文章主要介紹了Python3.7 讀取 mp3 音頻文件生成波形圖小編，本文通過實例代碼給大家介紹的非常詳細，具有一定的參考借鑒價值,需要的朋友可以參考下2019-11-11
Python for Informatics 第11章之正則表達式（四）
這篇文章主要介紹了Python for Informatics 第11章之正則表達式（四） 的相關資料,需要的朋友可以參考下2016-04-04
這篇文章主要介紹了pytorch中使用LSTM，可以在troch.nn模塊中找到LSTM類，文章圍繞主題展開詳細的內(nèi)容介紹，具有一定的參考價值，需要的朋友可以參考一下2022-07-07
這篇文章主要介紹了Python中的exec、eval使用實例,本文以簡潔的方式總結(jié)了Python中的exec、eval作用,并給出實例,需要的朋友可以參考下2014-09-09
python進行數(shù)據(jù)合并concat/merge
這篇文章主要介紹了python進行數(shù)據(jù)合并concat/merge，文章圍繞主題展開詳細的內(nèi)容介紹，具有一定的參考價值，感興趣的小伙伴可以參考一下2022-09-09
這篇文章主要介紹了python 畫圖 圖例自由定義方式，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2020-04-04
PIL.Image.open和cv2.imread的比較與相互轉(zhuǎn)換的方法
這篇文章主要介紹了PIL.Image.open和cv2.imread的比較與相互轉(zhuǎn)換的方法，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2020-06-06