Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
, F, e( ?' n2 p0 s- e( e
% p7 v8 j5 P5 Z3 C易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
' T4 X& {# C% N- r3 {! t7 g/ w' h移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
; {6 Z& U3 g' f' U$ L感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
  G( R% Q  ]0 d# l% K7 S同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
& X* W( A" f7 m7 a# N分别将时刻t处于易感状态和感染状态的个体所占比例记为，
7 J7 z1 E. G: m6 Q0 Ys(t)=S(t)/N
& ?2 L% `6 V$ c% \0 `8 V+ Ti(t)=I(t)/N
" }$ f+ B4 a9 s+ N/ s# i( U显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
2 y3 v& v7 {* E. c$ gdi/dt=βsi
) R8 e6 s8 N7 C. B. Y( j即
di/dt=βi(1-i)
6 h. Q8 l  \9 q, B: O7 U# }上式也成为Logistic增长方程式(Logistic growth equation),
6 [: U  o' c7 a8 D& i) Z  o方程的解和图像如图
" U3 E) x/ f* O) |9 H1 U
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448

5 O4 j" S6 r. ~- b7 W. A
'''
7 N% ~( L6 m6 {" H实验环境Python2.7.13,igraph包，cairo包，numpy包
, f4 p3 s5 c2 N, ?# ^% }7 [9 s'''
$ A% n% W: A. d' k- @# -*- coding:utf8 -*
+ \5 n. P% ~: b4 e9 kfrom igraph import *
import numpy as numpy
7 S+ y* T* i/ R( k) w' Pfrom  numpy import *
3 ^2 Y8 c8 p8 P, b, Vimport random
7 l4 ^: `. d) _5 e( o' X5 x
! A& b" F$ f) d  T9 Gdef len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
6 `. f) R7 w$ U, Q: |; H" Z, M
g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
3 {# {: T. h( Msummary(g)
7 `  H$ R( i) R. z9 t) W& l: gnodes_num=g.vcount()#统计图中的结点个数
+ T8 w4 o7 N; }7 w. E7 dnet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
! B) i, g7 ?/ |g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
. Q/ X1 c4 t4 f- g' X9 A8 znodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)

/ G8 Z6 \# T0 G$ Vinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
+ t& @' w1 d* {4 b) I$ \" G+ Tset_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
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
6 M0 K4 m( ~; G) _$ |infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

/ d5 f1 |$ @7 e  p" F+ I0 Hstop=False#感染过程结束的标记
temp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化
  R# m% h4 H. K7 U  G! u0 D0 N6 ^
while not stop:
# _5 J& Y5 i& q" u    i=0#记录让每个感染源都传播一次
    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)#获取本轮的感染源数量
        while i<temp_len:
, F0 A& J) @" s% O4 a( n7 B' L; [            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
# Y5 t2 e# N2 q( I7 i            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
) M$ s4 x& [  y( E9 G                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
- B1 f6 V, i5 B7 V. i            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
9 a0 w% A, [+ |8 ^* Z            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
: G1 Y  ^, Y% }+ n: ?                    nei_arr[t]=j+1
: }) M5 f* F, L8 [                    t=t+1
' j4 R  B/ U" d3 s9 z5 _            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:#存在需要被感染的节点
% a& X/ ~( M9 \* `4 S                t=0#让ran_infe_arr内每个感染源都被感染
0 W- N; }2 A  s8 S4 x- \                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
8 [# l2 x) {1 F                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
9 J, w$ v6 Z1 `! {' d                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
4 |5 F. p" V8 ~# n$ D                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
, y5 I, h- p( p, f2 p& W                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
, c+ ~/ k* C2 ~' ^. V$ j, Z2 F9 [                    t=t+1
% `3 C% ?; O/ E            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
4 S/ b. Z6 X2 p0 N" z7 n  n  |

视频演示bilibili传送门
效果图







& W' e3 k3 V" @$ a: t3 q3 Q# @/ U

' ]5 x9 E. y0 w7 U6 T

6 @0 d+ b% k- A
/ |8 ]  k( c' K( o6 v% [ ————————————————
2 o: \: x2 S3 M# v版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
; `: H, N; @* q) y! P原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
$ T& ]- n' U8 Y6 Z7 V& a8 v
: Y' o$ `% R1 I1 B. D4 H0 ~

尔雅

发表回复棒!
, t9 K5 t4 l- h8 G6 ~