Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
0 f* r0 Y, l6 c#SI疾病传播模型的原理
9 X; k1 j9 M" c& Q& [在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

! t  y4 V- o$ ]  y易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
* o; y0 e0 L7 r& U易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
, w9 a8 E- d- A/ R5 o移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
5 `4 o5 b0 v7 n5 D% rS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
0 B5 _8 {8 n. d  @' w感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
( n* f. s9 Y) W6 |3 `ds/dt = -β*S(t)I(t)/N
' ?% h: ]5 n; W7 B同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
9 x& Y7 L8 Z1 P& _! Hds/dt=-βsi
di/dt=βsi
即
di/dt=βi(1-i)
$ y" ~! L# h. F  U5 F/ S  z上式也成为Logistic增长方程式(Logistic growth equation),
/ Q9 {$ U  T/ X; V8 V8 w3 T方程的解和图像如图
  _( D+ z% [3 C  F
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448
+ v; l2 P/ }' [3 h; W$ b

  O/ n. i$ q: o7 U'''
) j/ {7 a5 |# t5 r实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
* b* g% p0 m8 `7 O7 w" h# -*- coding:utf8 -*
from igraph import *
; n: G% O! B1 N- _' Y; ^% Limport numpy as numpy
& r/ v# |7 N$ U: S) U$ h% tfrom  numpy import *
  {4 a! i  e+ T, i4 {import random
& s$ p* Y2 V; F* K6 Y+ Q2 b
def len_arr(infected_array,nodes_num):#获取感染数组长度
# u5 m# o# B) ^* U; i    len_value=0#初始化长度
- f: i2 ~) R0 P  s    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

6 o* k/ |9 _8 a' o% c. mg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
  R, K9 _9 l9 A; X" ^- hsummary(g)
7 V( P- a* y# |7 enodes_num=g.vcount()#统计图中的结点个数
3 M1 d2 ?- M/ q0 X* z# `net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
2 B- Q  @% f) {: u) \g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
; K9 n3 I  z) k* @8 _( y7 ]8 hprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
- ~. {- k4 f3 V! h; E$ C0 ]
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
8 D8 d$ j) h' \g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
7 b, e% p9 U- d% p
stop=False#感染过程结束的标记
temp_time=0#第几次感染
! Z. ~5 ]$ ^5 N$ Rtemp_len=0#本轮的感染源数量初始化
0 V+ b0 ^: a5 s% j. R
* G( i! ]; I0 S. }1 cwhile not stop:
    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
9 l4 v4 v, j4 }; a9 g& L        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
5 }6 `/ S! a3 Q$ o            for j in range(nodes_num):#遍历节点
  X8 D" G8 D- t5 \# N/ `( q                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
* X  @% o% k2 x9 O0 r            t=0
( g6 {1 y& {- m" v6 d* n            for j in range(nodes_num):
. b+ Z# A( J9 z8 e6 A: U% q# z                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
                    nei_arr[t]=j+1
0 k( a3 S9 T# o7 [. M2 ]8 B                    t=t+1
2 R$ z2 P2 k/ p) ]0 B            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:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
7 A# w2 M0 A# d" g. k6 A                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
( b  \; L0 [! z  J                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
9 ]& b- x: u" }( m4 {" h                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
5 S! j* J0 h0 P7 Y                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
1 [1 B: c2 H  z5 K3 `            i=i+1
* H: h% s3 K  b    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
0 K: b9 `+ R/ R0 I        stop=True
8 O, p7 |& P" }" ^
5 j! ^" J4 X: ^: o; A% m/ b
8 Y$ T  _( W7 o2 l# Y视频演示bilibili传送门
效果图

9 n* ~2 j, W; m5 ]


+ S2 n+ f7 n! ~. q, k4 T
, D! D, b6 C0 O8 X6 }: P( U
5 H5 o- \& x5 F, N
2 I2 N9 G" ]4 j' b
' ?- {* O- m& [) S( _' k; p- f
/ }. P  ^; l" \' I8 Q) r$ q% ]

; h3 Z0 q6 i2 [# w6 ]
3 x% W8 c2 k) i6 v, y8 a ————————————————
% o9 }( C- p/ ^1 U9 [0 C版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
% x4 p3 `/ ~$ o8 S8 q
7 }! ?) H  L  U! D( S: n$ X
4 o3 e4 ]! X4 c* B

尔雅

发表回复棒!
. x2 x5 H. k. T9 c) O5 {. G8 M2 R1 F