Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
+ m/ E' Y6 s3 T#SI疾病传播模型的原理
% K3 ^- K! q$ i  v在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

. w- [0 D) y; n; w易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
. x8 t# |# B/ K2 YSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
) e9 r$ z/ F! I, a6 a; k& h即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
9 f3 A( O0 D5 O感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
9 P2 _5 h2 M8 \ds/dt = -β*S(t)I(t)/N
  G4 y: K( p* L& T& e: W5 {/ f同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
8 |7 {3 l6 w! M! S分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
( s$ `/ L# h, n0 @s(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
0 x- u) T/ ]& V: P' u4 ]# pdi/dt=βsi
即
% c1 k3 z! Y+ ]% W  Ddi/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
1 _/ U: H+ r' F2 h  ]/ \0 W! F+ R提取码：z448
7 J# K9 d2 z; g6 S: O( @; N$ v( U/ m
# s1 |9 i4 r9 F' m0 |* v0 z
'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
8 k3 U$ {  \$ K- |, Z$ m' f'''
( f. x$ Q4 J# a3 B$ x  G# -*- coding:utf8 -*
0 T0 F( q8 Z/ g* @$ Gfrom igraph import *
1 [2 i6 z" {% C3 w6 Zimport numpy as numpy
/ ~8 _" J1 v% {from  numpy import *
import random
, s: {- Y% r8 ?6 a; R+ e' L' K# Z$ O
+ C% ~/ d% u. h$ L1 `4 S! Jdef len_arr(infected_array,nodes_num):#获取感染数组长度
4 [, p, J) L  `3 B. m- `    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
6 U2 i; p( A% }$ Unodes_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
  J6 G6 r% B& C, _- q8 W6 lnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
8 p9 O8 L" P' ?' X% _6 J                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
8 S! L3 F% g1 `5 Gprint(infected_array)
+ n1 ~3 b" E! a; R) _9 V# m
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
. s: R+ n8 d/ |/ N' z+ \8 e: ^set_time=2#传播次数(感染次数) 2次
7 O9 L5 t4 b0 R2 y! i3 |; `source_seed=1#感染源位置
4 n" u0 Z8 }# M2 T2 lnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
; d  w: Q8 t) p* Mg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
$ {/ ]+ K$ `' k& J, r4 ?* Tinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

stop=False#感染过程结束的标记
temp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

1 v0 [+ j% L: Q. S( vwhile not stop:
$ b; {( J, e9 J  \$ g# Q' N2 x    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
# U+ q7 I. q. N( J+ e        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
8 x1 ?3 e8 I5 c; b        while i<temp_len:
# U6 i% X" H, u* m5 A  B6 m" G            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
& j1 Y6 [: d1 T; V' S) m            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
% Y, l5 n/ x6 L, W- Y                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
" a( K, w7 B' e# q9 q" n                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
5 _8 o2 r  p# x. }            t=0
            for j in range(nodes_num):
( w+ W9 r" _* ?5 o5 I% ~0 N- _. H                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
. k9 ]- U6 z3 @: f% i3 B( f9 x                    nei_arr[t]=j+1
0 u) u& ]+ P$ a% h                    t=t+1
1 Z, s$ p3 v+ I+ ?( _            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
  n. v$ m* K' w$ K# @; [) G. Q                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
" U5 Z. d$ f3 T2 F: p9 q0 `            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
) Y# U3 U; w5 r) j9 q+ E7 w                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
' U+ w8 N; o+ O+ V1 z' b" M- ]                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
' e* x% F# ^  f/ [$ ]/ W9 U                    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)#绘制
. n: P! q6 a* C% [                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
! a, K" `. q0 \0 H        stop=True


# X# s8 E* X0 j3 y2 m5 C9 c视频演示bilibili传送门
. x1 j" ?4 t  L效果图

: j& ~4 e' k. E/ O( h/ Z



, }$ v# Q8 E! s

, g% d  u- T+ W, S5 m: K




. f1 X2 q2 N1 `) n" L ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
9 Q8 F1 n& A" H! ]- d原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

; y9 C9 E) a; n) [+ E, U  v

尔雅

发表回复棒!