用Python製作在地圖上類比瘟疫擴散的Gif圖

受傑森的《Almost Looks Like Work》啟發，我來展示一些病毒傳播模型。需要注意的是這個模型並不反映現實情況，因此不要誤以為是西非可怕的傳染病。相反，它更應該被看做是某種虛構的殭屍爆發現象。那麼，讓我們進入主題。

這就是SIR模型，其中字母S、I和R反映的是在殭屍疫情中，個體可能處於的不同狀態。

• S 代表易感群體，即健康個體中潛在的可能轉變的數量。
• I 代表染病群體，即殭屍數量。
• R 代表移除量，即因死亡而離開遊戲的殭屍數量，或者感染後又轉回人類的數量。但對與殭屍不存在治癒者，所以我們就不要自我愚弄了（如果要把SIR模型應用到流感傳染中，還是有治癒者的）。
• 至於β(beta)和γ(gamma):
• β(beta)表示疾病的傳染性程度，只要被咬就會感染。
• γ(gamma)表示從殭屍走向死亡的速率，取決於殭屍獵人的平均工作速率，當然，這不是一個完美的模型，請對我保持耐心。
• S′=?βIS告訴我們健康者變成殭屍的速率，S′是對時間的導數。
• I′=βIS?γI告訴我們感染者是如何增加的，以及行屍進入移除態速率（雙關語）。
• R′=γI只是加上(gamma I)，這一項在前面的等式中是負的。

上面的模型沒有考慮S/I/R的空間分布，下面來修正一下！

一種方法是把瑞典和北歐國家分割成網格，每個單元可以感染鄰近單元，描述如下：

其中對於單元，和是它周圍的四個單元。（不要因為對角單元而腦疲勞，我們需要我們的大腦不被吃掉）。

初始化一些東東。

import numpy as npimport mathimport matplotlib.pyplot as plt  %matplotlib inlinefrom matplotlib import rcParamsimport matplotlib.image as mpimgrcParams['font.family'] = 'serif'rcParams['font.size'] = 16rcParams['figure.figsize'] = 12, 8from PIL import Image

適當的beta和gamma值就能夠摧毀大半江山

beta = 0.010gamma = 1

還記得導數的定義嗎？當導數已知，假設Δt很小的情況下，經過重新整理，它可以用來近似預測函數的下一個取值，我們已經聲明過u′(t)。

初始化一些東東。

import numpy as npimport mathimport matplotlib.pyplot as plt  %matplotlib inlinefrom matplotlib import rcParamsimport matplotlib.image as mpimgrcParams['font.family'] = 'serif'rcParams['font.size'] = 16rcParams['figure.figsize'] = 12, 8from PIL import Image

適當的beta和gamma值就能夠摧毀大半江山

beta = 0.010gamma = 1

還記得導數的定義嗎？當導數已知，假設Δt很小的情況下，經過重新整理，它可以用來近似預測函數的下一個取值，我們已經聲明過u′(t)。

這種方法叫做歐拉法，代碼如下：

def euler_step(u, f, dt):  return u + dt * f(u)

我們需要函數f(u)。友好的numpy提供了簡潔的數組操作。我可能會在另一篇文章中回顧它，因為它們太強大了，需要更多的解釋，但現在這樣就能達到效果：

 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

匯入北歐國家的人口密度圖並進行下採樣，以便較快地得到結果

from PIL import Imageimg = 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 timedt = 1             # time incrementN = int(T/dt) + 1        # number of time-stepst = np.linspace(0.0, T, N)   # time discretization # initialize the array containing the solution for each time-stepu = np.empty((N, 3, S_0.shape[0], S_0.shape[1]))u[0][0] = S_0u[0][1] = I_0u[0][2] = R_0

我們需要自訂一個顏色表，這樣才能將感染矩陣顯示在地圖上。

import matplotlib.cm as cmtheCM = 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()

