Python实现简单的SI传播模型

杨利霞

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

易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
2 {1 ]' j5 f% L7 X. ^% \移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
* F! S, m3 V3 J* I$ \6 {SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
  S# w' o6 r: q, |9 ^在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
& r1 O8 e5 J- @$ F4 y感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
+ X$ j* L. t) M; H$ N1 {同时，感染个体的数量会以与易感个体相反的变化率增加，
7 G9 k! N' v6 v, ^$ ~' h% `) |" @& b0 nds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
0 B% j0 p/ A) E6 r5 u% Is(t)=S(t)/N
4 T; R8 r0 @$ P9 di(t)=I(t)/N
# t( k' y0 S" f8 _, ?显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
3 o' |  w* R6 ?+ kds/dt=-βsi
8 ^9 a' P% Z' U9 E; F( o% n) |di/dt=βsi
. _3 T1 t$ a( C+ T, U9 Y7 }4 l9 z即
$ P+ T2 q5 n' Y% k9 f2 ]di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图
- Q9 f: f. [7 q+ }1 N0 N
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448

$ g: l: t* Q' T* v# j
'''
3 H9 P" I- y; E% B$ R; h5 }( S& K实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
" X8 p. U  S! p! U/ ]5 B7 I# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
3 H, u7 j9 h( k1 p" l* r0 |0 @from  numpy import *
! B  T2 O3 T6 W: ^; Timport random
% Q* q$ _* D( Z6 l! i; @8 Z
& C7 i8 V6 s2 n5 g! mdef len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

+ `4 W3 }% c) Z( n' u; Fg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
nodes_num=g.vcount()#统计图中的结点个数
' [% ?$ t. o% @0 fnet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
% C3 g5 Z: a. m* C$ _g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
# a; y; J! V3 @, Y' C- `a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
% T; Q5 F6 Q1 p2 Jprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
7 g* w* [: }( f7 }4 S: [1 f+ aprint(infected_array)
9 z1 l% E7 L) k. F, C
- q8 f* P4 M0 _! vinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
0 M& t2 }, j: b' i( Nsource_seed=1#感染源位置
# ^' n& e. h; F  g& f1 bnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
7 P5 x. s* u: T* X5 qg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
& ~6 |+ Y/ t# a$ \- Z. z* R+ Hinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
  V+ U  H* c( W( Y# ]( i2 F# s% ~plot(g)#绘制

stop=False#感染过程结束的标记
temp_time=0#第几次感染
  o  p& w" }3 p# V) ^temp_len=0#本轮的感染源数量初始化
! y4 I' K6 B  J# T
# K! }% C0 I+ I$ N9 F8 i8 ]while not stop:
% W8 J0 [+ O! W4 k; p- K7 k    i=0#记录让每个感染源都传播一次
8 j, C2 ~; J' \5 W1 ^    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
$ K* x# Q( n! p/ e; Q        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
% K" o; e( e5 I) H6 g+ a        while i<temp_len:
7 I  ~! _$ l3 j) U            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
9 D/ A, G! H0 C3 I7 I  r) q+ \; M            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
9 f. _: k5 S& E& P* @            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
' ^3 I5 m& V% O  L                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
" F4 q8 b' j: |            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
! D! Q, B. z3 ?            t=0
# O+ i! _5 r/ \' R9 w5 G! P            for j in range(nodes_num):
& M9 Z  _/ @2 c7 R' ]" q                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
9 B1 o  \' z+ U  J- ]5 U                    nei_arr[t]=j+1
                    t=t+1
! [6 V8 D" y( f            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:#存在需要被感染的节点
0 X4 y; r% y2 v& x, w$ W' }  D/ x                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
/ E3 P; O8 n+ s2 Y& }% Q                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
; j1 `# I- B9 Y( t                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
8 W. I6 ?/ E3 v2 c9 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

: |' T. d  ]- g
2 E; D1 P1 c3 s% D8 x视频演示bilibili传送门
效果图


  W( h' m* T. _  E! N







6 l3 j3 y0 i& [

: V3 q' Q2 v9 G5 K2 U$ R3 I9 h3 z ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ [* [1 G- a9 n) i% I2 k原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
9 I: i6 b; Z0 S% i+ O6 Y
2 t6 i, r" |  E* P' }! W

尔雅

发表回复棒!