Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
% Q4 D2 y+ m( _4 z) L0 i, }4 S在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
6 x- x% J; T: n0 y0 [
易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
4 C' `" ]2 N, G" ]" V1 X* m移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
& r9 \1 Q0 R. T4 s! G5 Z4 o- A' u  ~6 iSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
2 h/ M2 {, T4 i) ?3 s4 y即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
5 r5 B, C9 q4 s" O. B4 Q1 a" u在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
$ x2 N5 H8 B8 g3 |/ G会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
2 b4 ^  H$ w& C- z+ M& n同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
% j3 W4 k6 P3 p; js(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
di/dt=βsi
即
di/dt=βi(1-i)
$ S8 T/ U* m, T1 k' v; {$ v上式也成为Logistic增长方程式(Logistic growth equation),
/ h5 a- L- x  K) \( E方程的解和图像如图
" \5 r7 \3 S% B/ G' j4 e
) }4 i, J( V3 `$ R  z' u代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
( l2 |$ G3 {4 n" |. ?+ r提取码：z448


'''
& U" x1 D0 k8 X9 g7 R3 i; T' S实验环境Python2.7.13,igraph包，cairo包，numpy包
9 O$ }; J! c  T'''
* [% B5 y: P+ r' U1 p& y# -*- coding:utf8 -*
, h/ {( A1 G2 {from igraph import *
' Z/ {' w3 y' _- g7 M. oimport numpy as numpy
from  numpy import *
  d0 N6 W4 n1 ~* himport random

def len_arr(infected_array,nodes_num):#获取感染数组长度
$ U5 J, N9 V6 P    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
1 {3 X! A  Y& Q7 a  n0 a6 }    return len_value
# I# I2 f( Q/ k2 D3 }
# o9 h& f0 R: ~4 C$ l) V+ Ug=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
2 Q- B% ?& v- U* r, Pnodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
  {, Q( W2 q) j) Cg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
' z# V% v+ V0 t, Q5 b" Ea=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
3 v) J3 }. e/ o+ E) o0 A+ ^nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
3 B/ a& A, k- c& l$ s) I                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
' `  D7 A" ?& Q$ L9 F$ J
2 x0 X% u3 u0 i4 X) ]infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
1 Q4 t/ J# l. z/ u# z9 o5 Tsource_seed=1#感染源位置
nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
3 F2 {0 m# C& E" E4 \# o! E* _nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
# c. V8 j  j+ A1 u3 i- Hg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

& u0 w+ m4 J" B8 _* @stop=False#感染过程结束的标记
4 g. w8 I' |$ ztemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

0 y! y" S3 {, {, Q" nwhile not stop:
$ \2 @( ]" O7 o9 E2 h& {    i=0#记录让每个感染源都传播一次
    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)#获取本轮的感染源数量
, w' m# H$ Q: e  `# O        while i<temp_len:
) s" ]# {$ h  F! l            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
/ T% {' |! A% h            nei_count=0#下一轮可以被感染到的节点数量
1 O% F. V" [& N/ r% g0 o            #生成下一轮可能被感染的节点的集合nei_arr
3 O: c: [4 j6 p1 J( j% u; K! x            for j in range(nodes_num):#遍历节点
" F6 B' K7 w* B  R) k. v                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
& c) l) @" a0 t7 A! l/ d            for j in range(nodes_num):
( k& r3 r  L* X& {                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
) {; Z4 B! v' `8 D0 w                    nei_arr[t]=j+1
                    t=t+1
( [  q5 B) l$ I% n1 _0 l            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
            if len(ran_infe_arr)>0:#存在需要被感染的节点
2 {' i  T: H8 i! W( p                t=0#让ran_infe_arr内每个感染源都被感染
6 a* F( f4 X: Y8 W3 j$ W                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
0 D$ ^2 p+ V7 |( [  m0 X7 G                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
! H3 V! y: J/ ^) B                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
" v2 F8 S4 b- x+ `# V% T                    plot(g)#绘制
                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
1 K# ~" h! A9 s1 M' o( Q        stop=True


  o+ m: ^' {+ `5 \! [视频演示bilibili传送门
效果图

7 {5 l7 T4 O1 R
' l; a% a5 i! d2 x: F; W8 }

+ M2 z: W" C) B' P- R) W
+ [4 R" \9 e4 g# a/ p! u$ I  ^





) g0 X( t% M1 h% L7 ^5 D  T
# u- P; h+ N8 O) m6 ^5 } ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

尔雅

发表回复棒!
# {, W: \. O$ W, `1 m5 C$ K