Python实现简单的SI传播模型

杨利霞

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

! v2 y1 T+ J' c5 G# s( _) a6 T易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
3 [. E  u, G4 f/ A8 ^SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
, }5 q' z8 c; v6 X; Q即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
" T6 t7 Q* V/ A$ A3 W4 P在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
6 m  x  }- O% ]( a3 g& ^3 w% O会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
1 |$ ^" O* r; ~9 T% h同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
# q; [4 M  M" @: K7 m显然有，
+ M2 Y7 V$ a  {( h' f; b. o6 \s(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
: y; m5 }: o0 n6 Odi/dt=βsi
" [9 R6 s& ~# H' L6 e& K. T4 G即
" g9 f2 \" O" U; odi/dt=βi(1-i)
; t1 n1 x! l  j  E+ S+ @上式也成为Logistic增长方程式(Logistic growth equation),
: J' ^9 Q7 l1 ~3 ?6 X- V: E  a& e# c方程的解和图像如图
1 I. L! l; k0 C" P7 d+ `  P1 o
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
$ f# f$ I% c: j$ z7 B+ o' K  X6 P提取码：z448
( E8 f7 j+ X! |' _7 k

'''
2 H& Y/ I1 T/ ]$ ^实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
# -*- coding:utf8 -*
) h$ |, Y! U) D/ T! L8 A+ @8 zfrom igraph import *
* r1 t: s* M6 \& Uimport numpy as numpy
, l' X# J3 X. z; J% A& Jfrom  numpy import *
import random
/ \. k: i' G: R  q  `
; V( J6 F" W: M# s9 b# c6 a/ e% `& mdef len_arr(infected_array,nodes_num):#获取感染数组长度
5 W' w& t# h5 I/ O9 H( X1 E& ~6 y4 @    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
+ ^% @1 z( f5 g    return len_value

g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
4 B( W; s; U* M( O' _0 ?summary(g)
nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
6 l  r2 ^- u0 q8 q: E2 ag.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
' h1 h9 L% E$ Y/ `; d1 [nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
+ c7 D) L7 {8 C: X# ~8 a& B1 A, S                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
2 ~# H! R! L6 V7 I$ b0 R8 sprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
6 E, ]' _) j$ t( Z# T7 {
! v# ^# `% t1 z% h5 I) Z% t' ginfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
5 _! y1 \/ n6 ~7 H1 Dnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
; k5 ~) f5 v! c% N$ n* Knodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
3 m' m) M8 i/ b# `0 T" Fnodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
& s+ U/ Y" _# X8 \1 ^  W1 sg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
3 p6 n' y1 D8 a# C$ |2 m( ~infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
% d( T1 m( ]: }' @
. {- [, z9 t( G; R% y( Sstop=False#感染过程结束的标记
temp_time=0#第几次感染
! n6 G! Q7 s: C# T9 K+ htemp_len=0#本轮的感染源数量初始化

' Y- [2 O. o- ?5 H( u  fwhile not stop:
    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
" `$ c6 M- Q5 q  w+ c        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
: C/ H7 [# F" _        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
$ q  ]0 g1 ]; _5 k0 c% k            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
$ @; m" ^5 _9 J) n) p6 [! b: ~9 Q% }  N. |            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
% C* X) N. c6 v8 ^+ `* f" |                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
) a; h5 X1 o. p! ]3 \5 a+ [            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
  j# U) g! \% ]$ S/ v            for j in range(nodes_num):
0 f2 h( {1 o1 g* h" Y5 h3 }                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
7 H3 [. |* `8 O  L7 T: N                    nei_arr[t]=j+1
) x$ b5 R2 `7 x' g) |3 V                    t=t+1
0 s* U2 ?2 D' c7 S+ a: R            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
& Y: b$ b( O: Z# L                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
, r( ?) n2 i8 Q& @/ [) Z* r; C            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
! t& ]! h6 w" L! D5 B" ]9 h- M                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
% U+ k8 s) r% |% o5 T  _                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
1 O' |( U8 ]+ s                    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"#将此次感染的节点集的所有节点颜色置为粉色
' }; ~" r0 A; H                    plot(g)#绘制
                    t=t+1
            i=i+1
/ S" X; U  }( h% c& S    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True

4 M; ]9 ~- R: H! d' Y
7 B* @, B' i% @0 O* ]. z2 B视频演示bilibili传送门
效果图
" v: w, q7 _* q5 A! T
4 R' g& Z1 {- c! x4 w1 t



. O8 C6 n& o5 O  T3 {
2 d; K& U* d" {6 O! C: e4 A


! C. O7 H+ e5 K+ j! f+ R. ^


————————————————
6 p3 |9 P6 z: A版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
" ~3 ?3 ?7 ]& }& W6 S/ |4 U原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
; k% I( \1 q6 L. F/ B0 v
4 W2 a% l3 y* i8 C
" j& t) B* l6 S6 E# T

尔雅

发表回复棒!