Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
1 F" C0 u; T  p5 [( v+ `; i" |在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

, s" [+ T8 I$ S4 h易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
. Q3 L. y  p9 _% Z易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
: E+ l5 T( g( e8 S7 V2 [SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
$ t- M) ]3 W; ]即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
+ G. l6 W. M# \" S. {) i在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
$ r5 M9 f  F6 O& F9 d9 a1 mS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
% z9 {# _: i2 ]2 N感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
! R% c- a! g  ~# Rds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
5 n% {+ t0 e) a& R1 R8 m, Q分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
* |/ n! X4 ?/ q3 F' r! ^/ u7 l显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
+ `" W$ {- t9 J. R/ z8 ~  g- w2 Lds/dt=-βsi
) l7 s8 i3 x! Bdi/dt=βsi
! v' M* g0 u/ ?1 o即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

+ e1 d' u9 Q  {+ ]7 O( J代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
1 Y& e! m- L4 A; ]. W1 f7 _5 u6 M提取码：z448


'''
% z, {4 Y' J- q: e  I" {实验环境Python2.7.13,igraph包，cairo包，numpy包
; v3 [6 k7 v2 p) L7 h- R& N'''
7 F  I- S! W# l3 y% ?& f# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
from  numpy import *
$ t, h# Z  u* C8 ^import random
5 y* e- M4 ?' G7 f
! S  X# Q- D6 j' @$ k5 S* d! odef len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
3 C, ?' O& I/ L0 z    return len_value

7 U# }) t! ]- {* ~% Y! w% _1 O: Cg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
4 t1 {0 Q9 D/ |# n2 l9 o, ^4 t. inodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
. ]9 y4 p2 [5 N) hg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
& o) ]; I( x$ H' ~& P& ]" e6 E5 Bnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
5 b2 t( b2 R4 X4 `3 Vprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
" q$ T0 G, u( v: D+ N( L3 `
9 V7 v1 |3 b9 P: C' l) Winfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
5 ], e, e6 A% Y% ~  y" Knodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
1 s8 `, `& W! r: V9 `infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
! L7 t9 O2 Q' c& q% }+ a0 x
stop=False#感染过程结束的标记
% }/ S. f. D$ Y* W- p# b. n2 d1 htemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化
, ^* e0 D, X3 V
while not stop:
    i=0#记录让每个感染源都传播一次
2 \8 ?0 O* E) b, L2 @    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)#获取本轮的感染源数量
9 S! ?+ d% ^4 b4 G3 X# |6 u. H& e2 d        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
( f" D3 l) ^( t, E( v' Q; w7 ~            nei_count=0#下一轮可以被感染到的节点数量
. c! @) ^5 ?. l  X% T/ f6 \: e            #生成下一轮可能被感染的节点的集合nei_arr
            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#下一轮可以被感染到的节点数量++
3 [2 l3 U+ k. k& H4 }            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
, l8 `+ ]% N) `: p5 H3 {            t=0
, K. ]. B, D( T" z            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
+ a/ E3 d) y4 ?- i- |, [                    nei_arr[t]=j+1
                    t=t+1
5 E3 B' c  J( ]  B" h$ X7 U            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
9 R- m5 B9 G& r4 o                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
            if len(ran_infe_arr)>0:#存在需要被感染的节点
$ H8 G5 F8 L* n5 a+ O                t=0#让ran_infe_arr内每个感染源都被感染
" R( D, H; M) h. s/ M6 |" E                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
- D5 t' H/ s( f) K! x% a5 w                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
/ ?8 b  Z* H3 l                    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]]#将此次感染节点放入总的感染节点数组中
                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
. n* I) b: T6 A! x' o9 {* w

视频演示bilibili传送门
( D$ y* i0 }* r$ a效果图

* z3 e5 u3 v0 @9 H% Y$ u6 R1 p
  U: }$ Z/ x- A# h3 `7 O
& k$ J. {: M* M! B

* I+ \# r( D% Y- W( `4 o8 G4 B, w






7 D, ?: d1 x5 o9 j/ u% e8 | ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
6 ?: M2 }) |% u5 Y5 V原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
% c0 N0 s5 B& \) t& g0 ]: r

* h8 N0 V- d5 T* y# s# g

尔雅

发表回复棒!
1 P; `$ o5 f6 F