Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
  r- q! u3 F7 H3 H; y9 }#SI疾病传播模型的原理
0 d9 D2 m/ H" [# m  K1 f- h- \  y& W2 D在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
9 X! R3 a" Q: s# z
易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
% J" f* {8 G- U即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
6 r$ d. U3 ^$ ~6 O' O2 i: h在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
6 C% M9 L  ?  i: @7 r5 v' b* dS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
  y5 U) T! p, H! u8 i6 d感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
2 P  \+ O' D; E- m+ Z: |  dds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
7 e* G8 v9 L; e$ `分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
5 R2 D: V: L5 w9 p. ms(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
di/dt=βsi
9 q3 f. ?+ n- K  {' b% U# D9 I; v( h* h即
$ o% O. d0 R, Odi/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
# F0 N- v. I; D% j4 K方程的解和图像如图

代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448
! Y/ R- I6 D- i( t7 E& Q' g. _

  C0 p; a0 ]" `/ u, A'''
$ w( R: g+ N' [4 I4 _6 n实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
# -*- coding:utf8 -*
from igraph import *
, e/ \6 w3 Q9 n; m# D2 l+ ?- Uimport numpy as numpy
+ @; q" v/ |& \+ }; @& J0 `from  numpy import *
5 {/ W9 }9 n0 z9 N# himport random
) [7 _* n: I! {
def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
) z! J$ ?' \7 p' z5 X" W$ l
, {% X; Z, A+ p  e( j8 rg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
& @/ A1 B$ R9 [2 p. K6 Rsummary(g)
7 `5 v* v& B; @nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
6 m$ t+ H1 Z" S6 z5 L                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
+ _) k! [' E+ C6 T9 S3 x/ gprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
; \2 P: n# P$ B/ @
1 b- v) V3 V6 [* w/ I+ oinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
  }0 J: u% w" L1 ]# ~4 {set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
8 |6 M. ]; g) w2 Z' Onodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
5 N5 {: B! v7 `4 z8 j) _8 M! h* Unodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
" F7 h; ]/ [( {  m. l% b6 @infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
9 |4 Z7 }- e9 Iplot(g)#绘制
  T7 ?' h! N' o
7 s, i: `* H: H: y1 N, V/ y0 [4 i# [stop=False#感染过程结束的标记
) ~, B" l# S: M2 O0 d' c  t6 `temp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
    i=0#记录让每个感染源都传播一次
) G7 N4 v2 V3 H, l4 k6 j    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)#获取本轮的感染源数量
" o2 T" T' {8 u' v( V% m+ G        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
* ]7 [% g8 J/ i! S  Q            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
. m" S6 j- G2 R9 o            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
7 q, M( g4 X) w& a6 Q                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
* Y$ s, g' ~2 P8 ?            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
. m% _7 v5 Z4 R0 ]- Y7 S; `, c            for j in range(nodes_num):
5 D* S) L+ i& Q$ y. v& M$ ~- H                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
                    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个数据生成一个对象
' Y, s9 e: h& |9 a            if len(ran_infe_arr)>0:#存在需要被感染的节点
$ @& [0 d% w" T" `% W8 i                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
; E1 `# `% o0 D/ R, n7 H! ~                    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
# s, M+ E0 \6 q& o  d0 ~9 a; \            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
' k! a1 J5 i' }        stop=True
+ Q  J# w. q, \- h/ l3 Z9 x* L+ s! U

视频演示bilibili传送门
6 }; s0 r/ P# o效果图

$ Y5 K- N: {/ z. Z# `  \
; s$ E+ t& r! b, h' U; W

) F* [1 f' ?7 O+ J1 a, @/ p, |1 n
  e' ?' ~. Y/ B% \

! i0 c! b- V, R" v, e: \7 v5 ^
1 L) Q( ~3 U* k5 p% r% d( `

9 U1 F/ n8 ^/ L7 T

————————————————
" G9 k* a+ u4 ]' c" H版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ ~/ H5 p& |4 l* ?; D. ]原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
+ w8 X) N8 ?8 N% T

. E/ m2 W$ d. Q# \1 h

尔雅

发表回复棒!
7 n8 z/ g( ]3 a% }: ^$ M4 C