数学建模社区-数学中国

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

作者: 杨利霞    时间: 2020-4-18 16:16
标题: Python实现简单的SI传播模型
Python实现简单的SI传播模型2 `. v* A- R! C$ x
#SI疾病传播模型的原理
; {& g( U, L6 M8 c! q# j9 s/ y$ |1 G在经典的传染病模型中,种群(Population)内N个个体的状态可分为如下几类0 |) h9 E8 D: Y$ v
5 b4 S/ p4 U) B/ p
易感状态(Susceptible)。一个个体在感染前是处于易感状态的,即该个体有可能被邻居个体感染。
2 i) M) ~7 }0 f$ {+ V, y: E易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。,即该个体还会以一定概率感染其邻居个体。- [; j& J* Z* ~' p) x0 U! f1 n
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态,当一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑改革提。3 f5 _; U! l8 t3 D6 B! s/ z
SI传播模型是最简单的疾病传播模型,模型中的所有个体都只可能处于两个状态中的一个6 }3 q( @! s0 g; N2 M: k/ O4 G
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
0 p9 [# U' X1 x9 n' P+ k在给定时刻t,令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目,显然有
# a" s# j. n& ?1 \9 _! _S(t)+I(t)恒等于N,这里,N是个体总数。随着时间t的增长,易感个体与感染个体的接触( [1 k/ j# Q7 R$ k
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β,疾病仅在
8 {8 a; c9 g+ i感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t,易感个体的比例为S(t)/N,感染个体的数量为I(t),一次,易感个体的数量将以如下变化率减少
' M' |( M3 M: M3 T- Hds/dt = -β*S(t)I(t)/N4 Y; ?% N% B' a$ s
同时,感染个体的数量会以与易感个体相反的变化率增加,
. F/ I! \6 v: Y3 ]+ P0 b/ \ds/dt = βS(t)*I(t)/N
# T8 Q1 t8 i# Z$ y( d3 u分别将时刻t处于易感状态和感染状态的个体所占比例记为,/ q% t! o* f) {# T9 x) n
s(t)=S(t)/N2 J9 G' R; z3 @# V' X$ ~* \+ t
i(t)=I(t)/N
) }) k! J$ d* Y" S; L6 o7 B显然有,
0 p% g3 ]4 d5 w- A, l2 u+ qs(t)+i(t)恒等于1,此时之前的公式可以记做+ z/ G4 L. h9 k$ }2 W
ds/dt=-βsi
! A0 r! |7 z; h+ x2 b. Ddi/dt=βsi
2 u9 R- f5 T" B$ h7 l6 h) ]3 O( E% `8 D% q' T4 r" T' Z
di/dt=βi(1-i)
  Q) q6 K2 j$ ^% D% Q) ?上式也成为Logistic增长方程式(Logistic growth equation),
' Q: a& [+ @# g4 [5 J0 \# L方程的解和图像如图
( R9 Q0 Q# F) \" ^* {0 P. s- |* P; ` 1.jpg
% O9 u9 v! X% v2 v$ b% J% n代码和相关文件以及环境链接:链接:https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
% O+ |' p! p+ G' H提取码:z448
. ~# E/ A. k" a/ C5 P
. Q: ^  R0 S5 s( I5 ?2 I  @* C4 X& `8 W7 J1 x0 Y, `6 K( d
'''
' O/ D! s+ L: e% M  c3 O) f4 k实验环境Python2.7.13,igraph包,cairo包,numpy包
' |/ }! E: N3 p5 }1 F4 E'''
, m0 R1 w" c5 a/ J% O9 D/ @# -*- coding:utf8 -*8 ]* X: _# q9 N4 R* Y& H
from igraph import ** a, ~) m& B! X& v  a6 o: [
import numpy as numpy
3 z' Y  O# q) `9 j- \- Q6 P# C: wfrom  numpy import *; L* C) f  w/ [' G! J) r
import random2 {: K: D' d$ R3 r; t( j) e6 s
; H+ {! m4 Q/ v5 g" @- G
def len_arr(infected_array,nodes_num):#获取感染数组长度
% ~6 S, V' u* k( q8 t    len_value=0#初始化长度
" p( n5 M& v# j1 W4 s( `# V& S. k    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
( Q2 P" ?6 \, |' m    return len_value
8 L" D) f6 B; g7 _( U# a3 G' i$ m" E* v7 C
g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
0 _: B0 ~/ B  U6 u/ n5 {summary(g)
( z) |1 W' s& H% L$ Snodes_num=g.vcount()#统计图中的结点个数
9 n8 E, o: f( u, ^net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
( K& S" T9 a. d1 G9 Yg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色! n% c' j+ n0 {! p; z. _$ U
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
. L, F4 U2 K0 P( p9 ~: Tnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建,用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重!!!)& \' ^( O2 z7 b& g5 Q- y
                                                    #第一列是节点编号,第二列是节点状态,感染状态用-2表示,第三列是节点感染的时间2 e+ t0 u$ l$ o  u, A! H7 K; c' z
print(nodes_state)/ |! _% d. t2 `' u. m
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
. N& \4 A' m) u; ?print(infected_array)$ e% |, M' D3 R/ Z/ E
. X4 H" k- P' {3 L; n
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染* v- G* h  y, `# w0 Q. ^
set_time=2#传播次数(感染次数) 2次
, w* L! c; x. O- }# s! L) E, ^7 Hsource_seed=1#感染源位置
7 U8 W8 [7 X3 u" Lnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
; x3 D: n- `7 H9 b8 qnodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态; I2 f) ~; `" x
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1/ C1 E' g4 K, m5 @. h
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
+ z4 ~, {7 u$ ^/ B7 Q9 W  c( X  Pinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
) W5 W% q7 F! tplot(g)#绘制" _% N' U8 M1 j
: Y" m. I' ~9 s# n' y
stop=False#感染过程结束的标记6 b5 l4 _3 B" u; l1 {  j, U
temp_time=0#第几次感染2 W% Q  a+ K' e
temp_len=0#本轮的感染源数量初始化
2 q( D* I) F. W8 O8 L) n1 i  t9 V) N/ e6 K3 ~5 _% j0 K9 w" ?
while not stop:
, }$ v: G* x$ {0 D: v    i=0#记录让每个感染源都传播一次
4 G9 i) I8 r) L    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行! }; q: }4 q' x$ P
        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量$ `( c+ W0 B* [
        while i<temp_len:
( G$ p9 b* {5 x4 b0 D' S( L            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间( D8 }" O% `0 ?- @& [; Z" b$ C1 m
            nei_count=0#下一轮可以被感染到的节点数量! {  R9 \& V8 c, [1 }9 ?
            #生成下一轮可能被感染的节点的集合nei_arr
! p9 d+ i4 y, k; \- u0 B! W8 t# ^            for j in range(nodes_num):#遍历节点
  P# L, ^/ p$ I/ q6 }                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染) j4 D, E) M3 ]' `- v' U
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++) M. M% X1 ~3 Y8 y4 R$ e. F" h
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染! \$ `# `! e+ B/ h# C
            t=0
' ]8 r/ m: s2 A: ?4 J            for j in range(nodes_num):
8 B, ?6 ]/ w: H* q6 [& ^                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:) [% B: v; ]  E( Y8 V' H; N
                    nei_arr[t]=j+17 k! s$ A1 }# [5 F1 L
                    t=t+1
2 q$ O( ~( R8 y$ i- S5 e            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
9 W( ?) E  D' V3 O- V0 y                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
) r. }2 I' _1 O! `( O/ K0 H            if len(ran_infe_arr)>0:#存在需要被感染的节点
/ V$ b/ H- u& B1 D" N                t=0#让ran_infe_arr内每个感染源都被感染
: \  k4 c) S# U0 N2 q" B                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
9 \/ a+ o) ^. _. l                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
( G+ N! h, z2 m7 a& q1 u) D+ k                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间9 i) G6 G! ^6 j
                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
1 b: v) q- j5 ?& z& C" }7 {1 c                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色8 I+ c; h) ~3 n+ p& K
                    plot(g)#绘制' P. ~# K$ R4 N6 F
                    t=t+1
& q/ i) |: b" ^7 v! R+ F3 K            i=i+1
( r" W) ~8 P  i7 i2 W. I" m  E    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
3 ]% R0 {9 u2 i8 u8 O        stop=True
$ w; X2 U: e2 ?1 Q5 t/ g/ L3 Q2 y9 c8 R) z) E

0 A3 S" V( z& p( {6 E$ ?7 L. G视频演示bilibili传送门* W+ l7 k* j2 v* U
效果图" r7 k6 ]/ ]8 k- Y, |
; A0 k+ |1 t5 f$ s4 M8 q
* n5 [- R3 j+ p3 d8 h  t: L! ~: i
2.jpg
4 _, h1 T/ y3 g4 e% w2 t4 P! c( i. [! g2 C  z5 v8 c  L: W4 F1 O3 h
3.png
. l7 _* i# d3 C9 h9 w( H/ M5 E0 d
4.jpg . r) a: D: q4 L+ _( {8 J

, F# Q' t9 f# k' `# r/ a) x 5.png
: }  s% n  h: C- S
: H, p5 x3 R, @) ^) s 6.png
  S; Z* j4 S- [# D3 L5 F# L- U, k0 |
7.png ————————————————$ M$ i3 ~- |+ I+ F5 s) d/ U
版权声明:本文为CSDN博主「eck_燃」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。) `9 b$ B8 }- X2 J
原文链接:https://blog.csdn.net/wdays83892469/article/details/80878862) c- m, l- y( V8 |+ `# T2 T1 J/ q

" N9 [" A1 H: s8 N8 ~% x
( F# U7 ]% S' h/ M$ H. F6 J
作者: 尔雅    时间: 2020-4-18 17:50
发表回复棒!
( Q1 n$ E* l  @/ x: ]




欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5