数学建模社区-数学中国

标题: Python实现简单的SI传播模型 [打印本页]

---

作者: 杨利霞    时间: 2020-4-18 16:16
标题: Python实现简单的SI传播模型
Python实现简单的SI传播模型
#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
/ \, O" W$ i4 ]
6 ~2 [$ e; C) `2 O3 H; p  x易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
9 Z' {: P3 X/ D9 B移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
* t* r0 w. ]- A3 ~: USI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
5 ~0 b, ]  D' _. A, L" q即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
, y# f/ m) W/ s, b  B5 a. w会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
' q2 c. U8 n! M- j' Y) _ds/dt = βS(t)*I(t)/N
" E- H" h/ O& g分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
/ |5 D4 i6 [$ i/ i% w$ x1 N/ Q显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
di/dt=βsi
即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图
! S) r! t0 `- Q, j# ]- q
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
/ H( w; e2 u1 r) C( V( g, d! ~提取码：z448


7 M# U; b  X. a/ l" J5 r: ]0 j'''
1 o9 o- F1 Q8 Q% y. {$ e3 c: H实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
# -*- coding:utf8 -*
from igraph import *
3 f1 O6 j: h4 k% g; D' o, Jimport numpy as numpy
0 W6 M. \- q/ U3 Efrom  numpy import *
- i4 k: C  v8 S: {! V2 T; E4 himport random

def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
! n+ f5 v# M4 p& Y. U/ h    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
3 U% b/ c% i' J% v0 E9 m
2 M1 ]" ~7 g! k( ug=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
9 L, m$ Q8 @; N. S" `summary(g)
nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
4 f/ [6 U& X, b+ `/ u4 r9 @a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
" N; @( J. L. }/ u; Inodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
2 C& t% D. h8 O; `                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
9 v9 P! I# L2 \9 S- C% Q, \) Jprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
2 d/ s3 F3 a6 V$ _$ X* Q6 P9 nprint(infected_array)
9 j5 o0 `$ z2 j
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
* J* V3 R8 S) O6 jnodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
5 V1 k8 @6 _; r; z4 |0 Pnodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

stop=False#感染过程结束的标记
% F' _/ ]- @0 F; {temp_time=0#第几次感染
, ]8 K" D8 g1 m$ v: H$ Stemp_len=0#本轮的感染源数量初始化
' N/ [% r# T* \0 L! V
while not stop:
    i=0#记录让每个感染源都传播一次
5 ~* ?; \, D, I6 D7 B5 Y    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
4 V4 s" q" D2 x        while i<temp_len:
! ?+ l8 F3 l9 V; Y            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
5 Z6 R+ t& r/ D! {, h% o- e            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
  B; f  g- T8 R- u( Q' E5 L                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
* L# `+ C, h  Y  g. Y9 P" [  D            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
0 A4 b+ e# t. {; R0 C            t=0
* G3 A: v  g! d$ B$ i8 r            for j in range(nodes_num):
; r. m! b( u8 R  k; r, {                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
8 a3 t% Z! {* E  Y4 l+ N; D, ^                    nei_arr[t]=j+1
                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
: n0 j* Q( P2 x, N                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
1 u. i; \: ~+ o# }4 M6 R            if len(ran_infe_arr)>0:#存在需要被感染的节点
6 r/ \( v7 h9 I7 S0 g6 f; f                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
- k" x; i) V5 C& V) E                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
: }: f3 v7 O: Z1 @8 h4 i                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
7 `9 ^0 j0 O7 @  R                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
5 A4 j0 m1 A" L% i7 h- M            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
! [. k4 ~% q7 }$ t( a% B5 y3 D
/ W* H* S) {) y4 I+ _
  U% E/ ~/ `( _. o# Y: v视频演示bilibili传送门
- V& [  F1 G9 Y' ^& O效果图
3 S/ o) Q0 x: M7 O  J4 \
9 R+ s0 m9 f. X, m2 j
* x) |) B: f) }+ o% G2 n$ w

6 B/ a& x6 Q! s7 M) H1 v
+ T# q- z: Z4 v7 L3 J


, l6 j. Z: F+ N7 L$ q& O
& J! ?! J3 P% `2 w$ @8 s. R: u
2 f, Y, U# U/ u" i3 v
1 H$ l( b' y  d7 i' k0 C' \
, ?: m( x) r. Q' q ————————————————
7 B; Y- j, b+ h4 w6 T版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862


7 n2 C% H) _% J  m; k  L4 J$ c

---

作者: 尔雅    时间: 2020-4-18 17:50
发表回复棒!

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5