Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
; W( S+ k! p7 G在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

4 O& ?) t& _2 t& p7 {  e9 M易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
& _* R4 G1 x  m, t) S易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
" i7 T0 }/ L' ]" X% A* L/ `即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
- O1 A+ s, q6 Y; MS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
: [9 M+ U8 N. K3 Sds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
3 H# F2 c) n+ Z1 I5 Eds/dt = βS(t)*I(t)/N
6 r# b% T5 r0 Z3 \分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
3 J( a8 U. g3 j# s$ J; N, Nds/dt=-βsi
di/dt=βsi
% L% W5 w8 }$ Q2 n$ j即
) }0 x% Q  R8 m) t* j( Sdi/dt=βi(1-i)
$ |5 i2 u: g& c6 L: H8 c9 ?上式也成为Logistic增长方程式(Logistic growth equation),
( S/ L, r6 ~" l. D+ M$ k: I方程的解和图像如图

) g. H+ A6 I% @5 Q" n9 _. Y代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
& Q3 Q3 J( H; h0 H2 o提取码：z448
% E2 Z$ \) F/ Q

'''
: |% h$ E" b  m实验环境Python2.7.13,igraph包，cairo包，numpy包
0 F' V2 d( r& U" Q'''
# -*- coding:utf8 -*
8 U- C* `  t3 F3 e* `from igraph import *
import numpy as numpy
- w1 P2 p7 C, H4 a$ }from  numpy import *
% Y$ z( Y- Z& T& Q  `4 v6 v, vimport random
# H* ]- j- M+ G+ w- O+ E
def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
: C1 _! h1 q. D    return len_value

$ m; A4 f8 u& ?9 ~  `* b  C" I8 a( mg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
8 P; c9 V; u7 S# Dsummary(g)
7 @0 M) u7 d, O1 B; z& xnodes_num=g.vcount()#统计图中的结点个数
) T. I0 X+ x- S6 fnet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
6 E5 R7 `3 h% Y4 Ua=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
: A( _. C  C8 l# |7 q/ K" _nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
1 E# y. r) L, U4 d* a, Zprint(infected_array)

infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
4 e( n) c4 j- x, X. h7 _set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
- {* c6 m' l+ dnodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
3 }: \% P+ I4 V, n  B5 n' |4 ~infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
1 |8 }& D% ]$ q9 iplot(g)#绘制
2 Z, W5 k: @$ X/ i6 V6 w
stop=False#感染过程结束的标记
temp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
! l3 l/ S. Q; c0 B5 X0 u: A& _        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
% M& y2 b1 b7 c: Z0 g        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
$ \2 M# |; c; v7 Y7 D            nei_count=0#下一轮可以被感染到的节点数量
( I$ R0 Y7 K2 F# I7 E3 U5 [            #生成下一轮可能被感染的节点的集合nei_arr
" z! d1 n3 h; S( z4 ?            for j in range(nodes_num):#遍历节点
2 W5 i5 G: N8 k' O! W7 N/ n                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
4 K% a( d7 u0 A8 z$ \8 B4 @% k# h                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
6 W# r! M- ~9 B1 A' N0 ~            t=0
$ r/ _3 O. {% g2 V) Y3 G            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
  p* m- e* }/ Y% g                    nei_arr[t]=j+1
& C) o" H- W/ r: P                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
7 P- ~  d. G2 c' }* Y                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
            if len(ran_infe_arr)>0:#存在需要被感染的节点
7 u/ l5 c+ a/ c7 j0 I- L                t=0#让ran_infe_arr内每个感染源都被感染
( Z4 t: P* U* R# J6 {  B3 p                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
( E2 r9 C9 Z6 U' ?                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
5 }9 j; r) D# Z                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
3 g! m5 z) q% ]/ F0 g( u; u7 |                    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"#将此次感染的节点集的所有节点颜色置为粉色
) f& O( R  e5 K$ ~                    plot(g)#绘制
                    t=t+1
            i=i+1
. M3 N2 L/ u& Q4 g& F# g" `    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
& S, U! Z9 k  F3 i" t        stop=True
# ~9 p' O2 a* p$ C

视频演示bilibili传送门
/ F$ m( q; K4 S* H- [7 }9 V6 u效果图
: L0 M. v3 C# O' |) L



3 q0 a! N) f) @1 h9 }


8 H; ?' J5 ~2 t3 F! l: n8 _
: k+ q4 t- T4 |9 ^

6 J8 L( t- _, [+ `

————————————————
) S. a, t5 R5 M! f& ?5 I版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

尔雅

发表回复棒!
( O/ s8 L5 t, D! B( D1 ~1 M