Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
! X6 Y/ t3 \$ Z5 o' q#SI疾病传播模型的原理
1 `+ @% _; m' f在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
" G: n1 W" V$ F& q9 H3 M
易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
: H) l$ ]& i1 q( C6 \即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
- K$ K- v0 b9 x- q# J6 \' L  D% J4 h" {2 E会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
: Z' k3 k4 O* g. m3 v感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
( n- Q! m1 h6 n. k同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
- @6 W# P) L% _) k, _, ^% h& J分别将时刻t处于易感状态和感染状态的个体所占比例记为，
% y6 S2 m& [  p" fs(t)=S(t)/N
i(t)=I(t)/N
8 s! G; c( g6 R' u- w显然有，
% x$ `1 W: b7 z/ v* ]: H3 ds(t)+i(t)恒等于1，此时之前的公式可以记做
9 b( n% J9 }8 E+ Y0 F3 Xds/dt=-βsi
6 H+ O$ ~- \7 edi/dt=βsi
$ {- [. b) t, @* E& l; p即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
2 Q0 `$ b9 E( r7 b% v2 F/ |, z方程的解和图像如图
+ @# V, }# G: e+ K, a5 \  W' d
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
( j* \2 n, p. P- t- j! j/ }提取码：z448

2 z5 `3 Y( x5 \5 n) c, M7 f5 a
'''
3 ]" x. s9 ~$ d8 w' b7 L) j9 ]' i% G实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
# -*- coding:utf8 -*
# d) D. U, q" G7 w2 Afrom igraph import *
3 K3 \0 G9 M2 x) y0 m6 b1 m- gimport numpy as numpy
- E, T7 t5 W* T  efrom  numpy import *
8 `+ U7 F  Z+ vimport random
- g3 j% l7 V. h8 `4 P1 [) a
- ^) j% l) @% o" e$ sdef len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
% J3 Q! @# N0 U6 p  b) H    return len_value
( A+ o6 X- j# \; ~
& W1 I1 \: h5 F' M7 mg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
- z8 W, X9 Q" da=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
( p* c5 S; w! l! I) Y+ s9 m  u$ B                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
0 k- S* s% L# R% f1 f9 ?print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)

infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
8 p: Z2 s; d7 {7 |source_seed=1#感染源位置
3 _9 f! U# Y! K* l" _nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
3 S6 [- y2 e, `nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
  K% u. n& Q3 w( A+ lnodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
2 n4 E0 \' {8 O( hg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
( E. n& v" O4 z. f+ |; P  Kinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

stop=False#感染过程结束的标记
temp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
$ P( \) @0 I' v. P/ ^' Z) C" B    i=0#记录让每个感染源都传播一次
4 g; i, k& E8 q1 V, O) o1 X- G7 w! d    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)#获取本轮的感染源数量
        while i<temp_len:
( V- i5 c: _& l4 D9 d9 O% A            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
" ?1 i- ~, d# m            nei_count=0#下一轮可以被感染到的节点数量
. `, V$ ]- M5 S- P            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
9 k# Y* o8 m+ Q: L' [6 p6 b                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
: s$ G& h) [0 e5 g, e  K                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
3 ?2 B9 _$ Y" x& ?) l- k            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
. H3 u! [  p) |' z) O) [                    t=t+1
  Q4 g; o  d6 f$ }, }: u            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
# k7 d6 J; ^/ n% v! z& ?; u# G# G                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
  j4 u1 V" G1 `9 ?; I2 X            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
5 s; Y9 B1 y+ b% W1 X                    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"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
4 ~! W) x$ L; T6 z                    t=t+1
            i=i+1
- Z! k7 L7 l1 c. _( }- D    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True


3 y$ w4 \. s" ~  m$ W6 b: C- O视频演示bilibili传送门
效果图
7 o3 v" g$ o1 i2 d
1 r) @3 @- Z2 L( {/ m) ?0 P0 T
7 l* U- I9 w3 {$ T" r

. M* {$ q; a1 ?1 }5 I
! y- X8 e2 X9 s. [3 w0 a5 G
0 V+ Q  y/ @" y; N* e. X' C
( W/ r- D% U8 {; O, o% Z. A' V9 s
  K7 U0 e! r4 b' @

* z0 i8 G; U7 ?, x( {2 G9 M3 u: |5 J
5 h5 p! {' E/ A( `6 Z) Q3 g
————————————————
' M9 ]% Z# x. q( [/ E版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
2 |/ D) n; t* v% x( q# C

6 ^" A* p, m; E# k

尔雅

发表回复棒!
: c/ Q! ], j! D