Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
% S, q5 E( x$ s# Y0 D! x# j2 J#SI疾病传播模型的原理
: \$ s' |! N- y, e: Z0 \! }$ `& P; l在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
0 f' U9 U) t. }; f2 N
2 G  p8 x/ a' _( Q易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
2 e2 W6 v+ ?8 i! U0 G易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
! W) Z- a. K. B: E移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
- M' K/ x, n5 |4 y; k/ z- Q在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
! \7 V( e" y8 SS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
/ A1 S- S9 ?" S/ i+ d. G会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
+ D6 O' B* Y' t3 l! ]2 `同时，感染个体的数量会以与易感个体相反的变化率增加，
6 z. }! y( w& z5 h) W% [: |ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
8 p5 u0 T5 L; x# rds/dt=-βsi
di/dt=βsi
& A+ F  \. `. `7 W* g3 e即
+ x: v" j& V8 c2 K4 Bdi/dt=βi(1-i)
+ l$ P( W0 h! @3 p  {上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图
% Y1 }  v, u* J" u$ Z
' T9 O! a- _4 K5 F+ C代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
- s% I' O6 Y) c( g: k8 A6 v提取码：z448
# p1 m# o. C7 Z" [  O" r

'''
  P/ F$ `* s4 ]实验环境Python2.7.13,igraph包，cairo包，numpy包
$ P9 {# |. q( p3 B0 v! P8 X'''
# -*- coding:utf8 -*
/ K$ ^7 I7 W1 Z# O; qfrom igraph import *
6 x- f) X4 w: ?import numpy as numpy
- h- n2 Q( `, i1 {) Z& @from  numpy import *
" l) K% L$ z+ F" R0 {1 dimport random
- H; v  V/ e* u! B( ?
def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
5 o5 F  c1 q1 |' ?( z5 |
g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
, p/ D: s" v5 k0 n! \nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
- W$ \' f% ~- \+ R0 e. fg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
" K$ I$ e7 Z: i7 M; A) h                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
, ^  b5 ?' |8 t$ G0 uprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
0 T0 y; A- ^2 b3 K; u; E
$ |& `; c5 D( Q( _infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
+ G; C' O4 C( g' B5 X1 Qset_time=2#传播次数(感染次数) 2次
8 ], n( F1 @% j5 S' lsource_seed=1#感染源位置
* Z' h( D4 U" `6 w$ Y6 {/ gnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
, l$ [* J8 o% M% snodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
, l4 E  [' Z! B8 A0 l  vg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
8 F8 J2 R+ b: f6 B! G
stop=False#感染过程结束的标记
temp_time=0#第几次感染
1 q# ]3 Q2 E( p4 L) ~3 j, u, ~temp_len=0#本轮的感染源数量初始化
- o: i" r$ w' Y2 ?7 g: l  f
  ]( Z' g# l& ]( iwhile not stop:
' w% n) E  L6 y' ~5 m5 ]    i=0#记录让每个感染源都传播一次
; X. F6 {0 a* f9 }5 X; R; \    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
4 D: z: M8 H( ]        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
% s1 \% E1 N2 d- o9 I            nei_count=0#下一轮可以被感染到的节点数量
* W0 d! B$ S; w, c9 ~! M            #生成下一轮可能被感染的节点的集合nei_arr
9 N. O* H5 S3 J! Y. D* u            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
5 @3 ^1 l: Q, L( {& S* a            t=0
            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
                    nei_arr[t]=j+1
                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
7 p- o, v- S) p: ^" s  Q, p                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
+ H% u8 h/ }+ {$ W! a6 i4 l            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
! ?  k, c9 [) q. G; O                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
) u, z) {. k6 n0 P                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
$ S& Q+ h6 U7 w( |                    plot(g)#绘制
8 H; X; I  f6 D2 J9 l, e. G3 r                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True


; q9 B1 M! P* W9 W: f# H6 B- O视频演示bilibili传送门
1 l- u4 P1 |/ k6 D效果图
& K( m6 y  |( F) w3 U  q
% F% }5 D7 Q8 h7 {

' W. o; O' _: ~% {7 ^0 x




3 u; q$ _8 P8 V
1 a' Y; M4 c& M; [0 J5 N6 i


7 b; \$ |; e: ^+ \* P) S; u) V ————————————————
7 J. N4 f4 w) Y2 M/ z1 X版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
9 G* O- E6 P0 u原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
1 c" H$ L% }+ G) g- X* }1 a+ L

尔雅

发表回复棒!
& x( h8 E3 R, Q" z+ \9 p