Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
3 x# z* w7 ]( M, O#SI疾病传播模型的原理
+ {0 X8 A: t! }在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
; F. G' o" \* V# C0 O
% d, m# R" w" U, q8 b) Y易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
# M. K- W% p  M# B3 R- d. g! `$ b. t即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
6 p8 v$ k5 R! f+ G" A% G( t. SS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
! r2 M6 B1 q6 o$ l" @会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
* ~3 y2 }) \9 u5 R% Qds/dt = -β*S(t)I(t)/N
" I5 n  F; E8 B/ m% B同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
. t% f5 L) v* Q/ F/ W分别将时刻t处于易感状态和感染状态的个体所占比例记为，
/ T! u( q1 z/ e4 U& `) e( U4 e+ A# e% es(t)=S(t)/N
& h1 G: ]* |) ?+ ?; e. Ei(t)=I(t)/N
显然有，
" }0 ]. O# k: o" h2 p* E* L  e2 \& Os(t)+i(t)恒等于1，此时之前的公式可以记做
+ v( o4 s0 [! n1 Z1 y: F, I8 Yds/dt=-βsi
di/dt=βsi
即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

' q8 L' x, X4 T8 o: ^8 ^1 D代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448


% b- \9 n/ j* y'''
  B6 |8 A# r  o实验环境Python2.7.13,igraph包，cairo包，numpy包
1 O' M; X6 u9 r; |" ~, n/ r+ A'''
# -*- coding:utf8 -*
from igraph import *
" K8 n$ R6 Y% ~* Dimport numpy as numpy
from  numpy import *
import random

def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
4 @3 @( w$ S5 {4 ~7 g    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
/ M- }% G2 n1 x; _! r4 l    return len_value
0 w( z3 j6 ~$ R6 T
% K9 O! Y& c8 {4 a/ U- W; f$ g) pg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
2 F6 t6 r5 F4 y  g" `summary(g)
) g9 ^1 o, d" \& K, `nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
/ T6 ~7 |4 y' n- l6 o" Ea=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
) n7 Y' Z$ K  `/ F' i9 nprint(nodes_state)
) K0 J; I9 A9 i( ]infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
1 N$ w/ Z5 h# M5 K, eprint(infected_array)

$ V" v& L1 @9 s! [infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
7 Q0 A/ I9 {8 R3 U/ p9 Bset_time=2#传播次数(感染次数) 2次
8 A& w2 o9 I2 W% Ysource_seed=1#感染源位置
/ Q3 A) I5 O1 P/ a2 V, Znodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
* b* D+ p. B) I0 t5 f# inodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
2 }/ i' c8 t# \9 ?, J* r* enodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
% n, t% u/ d$ ~, s" \infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

! R6 W3 N# O2 W7 bstop=False#感染过程结束的标记
temp_time=0#第几次感染
6 o, E) ?6 v' m3 B$ ]6 W+ Jtemp_len=0#本轮的感染源数量初始化

while not stop:
    i=0#记录让每个感染源都传播一次
' N  k9 ]  _2 Q) P- \    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)#获取本轮的感染源数量
: J' \4 p! n+ f) e% E        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
( Z! i2 R* b2 b; W; U/ O5 ^, a* P                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
! H. C/ @! P3 `. E6 E: r            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
& R6 k: @7 P+ D3 z3 l                    nei_arr[t]=j+1
                    t=t+1
            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:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
7 H& i. O  P/ l% a* |% e: {                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
/ O" X* w+ w8 J8 V/ f                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
. o$ v# f+ |- C7 e9 \3 D                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
- |9 M3 u& N) [, A                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
- g" ?& z8 r1 t5 ?" R, h                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
- z1 S5 z8 L1 ^) R1 I2 d/ ]/ T                    t=t+1
9 f% q% @4 K. x/ N3 A1 Y+ Q" `8 c            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
7 x  q7 M. f+ f/ {

! y/ M5 A" e, ~# s; `$ z/ i视频演示bilibili传送门
效果图
% O" x3 }& d  L7 Q/ `; s) n
2 X- u6 b" ~! u' S( Z


# _  ~. J) J0 x


$ ?! j; l7 c9 R9 m$ c/ B

% g6 @, B4 g) e! r" i1 e


, m  K/ R5 Y5 U3 Q4 M ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
" f8 r$ v$ I+ c0 [  I8 q/ S5 Q
, s7 F. u7 L* E/ B2 \' f* c

尔雅

发表回复棒!
) h5 g- s- C6 j- g$ w