Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
: z8 k$ Z8 `# G# }8 ^9 W在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
+ L4 [' e2 Y) \8 A" d. h& I
; Q. G7 B: X- N/ A9 p" ~0 [2 z易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
' w% h, u' v. t- B% ~移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
6 ~3 S4 y- d% HSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
) C; X* {  w7 L" }即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
* O! T1 J+ P, r在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
2 @" d$ I0 y1 s0 d2 I/ H8 k$ FS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
' o) ^- n; d9 a& [& ]- j  A+ m" Cds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
$ x$ D, Z# N# |s(t)=S(t)/N
i(t)=I(t)/N
8 |1 l0 A0 |* e) }显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
* B1 h& ?) d( \2 }' i8 u$ u0 lds/dt=-βsi
% G  c9 J* W5 i5 z: {% Udi/dt=βsi
3 h; s. W9 }8 }/ \7 `0 y即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
2 B/ P1 h5 }+ J( h+ y( s4 q( }2 _% O方程的解和图像如图
" f: B5 m2 v+ A. P2 T
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448


% O5 D" m9 c$ \) F'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
! h% y0 d2 l& h; P$ _# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
from  numpy import *
import random
. l4 M7 Z$ z+ A0 P( U
) ?7 k0 I. B% b0 C5 sdef len_arr(infected_array,nodes_num):#获取感染数组长度
6 G, c4 [8 ^& z9 T+ D  z5 y    len_value=0#初始化长度
* k- w' u) E) `9 i    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
; t8 M" }$ X/ v2 M  r) d
: e0 ^$ `; W( @  @) {" M2 h, ag=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
7 i0 Z. b9 f5 B1 m; ~3 r1 ~a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
5 y" w' _8 i( A( r4 V& \! c5 |$ ]                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
" T3 H3 D9 f7 F/ N* x% N$ T7 Sprint(infected_array)

& g9 e/ G" u: Ninfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
* e- O( ~4 ^& _# u4 r% Wnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
# ]  [5 d2 K0 F. vnodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
- Z. p) X8 C5 E
stop=False#感染过程结束的标记
! I. w( Z$ H" Q) \1 h' Ctemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化
" `# w0 Y7 V) [( ^, y3 l$ {: `
while not stop:
9 m6 s9 m  V2 g' \    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
, J/ G/ y5 G0 i7 m& j        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
; M: A  ]: o' d1 ]        while i<temp_len:
8 {6 [1 t+ N! W+ e; O            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
, i" F6 E! c' l4 e# v            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
: y  I9 b) R) z8 n  e$ A, g2 l                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
2 [2 q( }. A+ T6 Y0 r) a            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
# |" z) `; m# S) Y* V: D0 M1 o! ~                    nei_arr[t]=j+1
9 _& |9 Y. f, I1 G2 \* E7 b                    t=t+1
9 i9 x8 V9 F3 v" m5 m            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内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
1 }8 d' D$ ]6 }& f) Z8 @                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    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]]#将此次感染节点放入总的感染节点数组中
3 b; o/ o% x( i9 L' h- ?" K9 B                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
- n  l7 P, z# ]+ d- G2 Q                    t=t+1
            i=i+1
% d8 I4 b9 E' A% Q  a: @' R1 x( Q    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
. }8 n; z, a! h; ]! ~% w7 C8 }        stop=True

, X/ k7 y5 E( V2 |8 i' [
视频演示bilibili传送门
效果图

: R, ^4 c1 u& p; C* i. }$ ?

- n, L0 ?+ m1 ?/ T

% j% q9 m8 [, f, N1 w# E% B

; |2 K+ H% T2 S

. r( T9 s1 q6 S* e' r

) O2 H0 e/ G* O# |: U* _2 ?* ^
5 d: D4 L4 _1 w/ W ————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
( c5 O# J. I/ G  U$ y) W% ^原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862


+ Q  H1 `% k' K- P

尔雅

发表回复棒!
% v8 a/ U0 m9 S! h7 ~