用Python制作在地图上模拟瘟疫扩散的Gif图
作者:Max Berggren 发布时间:2022-09-26 17:20:22
受杰森的《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()
猜你喜欢
- 下载了一个小型的记帐软件,发现这个软件数据库用的是access,很想看看它的数据库结构怎样,结果人家加密了。access的解密小case了,
- 【原文地址】 Recipe: Deploying a SQL Database to a Remote Hosting Environmen
- 有时候很多朋友访问自己的asp或php等程序页面时出现一些错误,就是页面无法显示也没有详细的错误信息,就算iis中开启了显示详细的错误也不能
- Python中生成器和迭代器的区别(代码在Python3.5下测试):Num01–>迭代器定义:对于list、string、tuple
- 前言1.装饰器本质是一个语法糖,是对被装饰方法或类进行的功能扩充,是一种面向切面的实现方法2.装饰器可以分成方法装饰器和类装饰器,他们的区别
- 在用JS编写动画的时候,经常用会到布局转换,即在运动前将相对定位转为绝对定位,然后执行动画函数。下面给大家分享一个运用原生JS实现的布局转换
- 1、需求分析首先我们打开地图搜索“鸿星尔克”:复制该链接到浏览器,发现这是一个json格式的数据集。
- 描述给ChatGPT的描述内容:python在桌面上显示动态的文字,不要显示窗口边框。窗口背景和标签背景都是透明的,但标签内的文字是有颜色。
- 大家好,今天才发现很多学习Flask的小伙伴都有这么一个问题,清理缓存好麻烦啊,今天就教大家怎么解决。大家在使用Flask静态文件的时候,每
- python查找多层嵌套字典的值def find_dic(item, key): if isinstance(it
- 在进行Web的交互设计中,颜色信息的传达也是不可或缺的一部分。我们常会发现许多“灰色”的应用,他们的出现总是不动声色而又恰如其分,维持了整个
- 作者:Scott Gerber原标题:Mobile App Development: 10 Tips for Small Business
- 下载安装Anaconda下载地址如下,根据所需版本下载安装过程暂略(下次在安装时添加)下载安装Pycharm下载安装Pycharm,下载对应
- 一、什么要备份数据库 ?在现实IT世界里,我们使用的服务器硬件可能因为使用时间过长,而发生故障;Windows系列服务器有可能蓝屏或者感染病
- 一、安装步骤 1.官网下载安装包2.安装一路next即可,安装位置可改到D盘3.添加环境变量将如上路径添加到系统path,不会的参
- 一维插值插值不同于拟合。插值函数经过样本点,拟合函数一般基于最小二乘法尽量靠近所有样本点穿过。常见插值方法有拉格朗日插值法、分段插值法、样条
- 1 简介费老师我在几年前写过的一篇文章(https://www.jb51.net/article/243348.htm)中,介绍过tqdm这
- Linux 自动备份oracle数据:曾经有个同事,来回操作开发和生产的数据库,结果误删了生产的数据库,那种心情我想不是一般人能理解的,虽然
- 在我们日常接触到的Python中,狭义的缺失值一般指DataFrame中的NaN。广义的话,可以分为三种。缺失值:在Pandas中的缺失值有
- 数据加载、存储与文件格式pandas提供了一些用于将表格型数据读取为DataFrame对象的函数。其中read_csv和read_talbe