用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 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
适当的beta和gamma值就能够摧毁大半*
beta = 0.010 gamma = 1
还记得导数的定义么?当导数已知,假设δt很小的情况下,经过重新整理,它可以用来近似预测函数的下一个取值,我们已经声明过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
适当的beta和gamma值就能够摧毁大半*
beta = 0.010 gamma = 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 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()
上一篇: 西兰花加工怎么做好吃又营养
下一篇: 板栗是干果的一种吗