Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
- \) U! x; p; Q6 `7 p4 Q. u在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
) O( ?# F# K7 @" z; _即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
4 V, M# O( ]& J& mS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
( s% u1 X9 @7 `5 a1 n$ m* M# _感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
! i( W7 X7 J, ~- |5 [: j! b% H  @( zds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
0 w$ q8 M- p  C+ p# M0 \s(t)=S(t)/N
" M; w7 y+ V1 x4 Li(t)=I(t)/N
显然有，
) T" y* h  u$ a9 s! f: Js(t)+i(t)恒等于1，此时之前的公式可以记做
  D5 l1 y# v7 {4 [& b3 `ds/dt=-βsi
di/dt=βsi
# j+ t' e* y6 ]8 m( B5 }; I即
1 ~' Z% N: o* N: m  qdi/dt=βi(1-i)
% O. K: x% i: L2 ?6 Z. p9 {( }2 B上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图
  t, G' Q$ K- e- O
代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
, w; Y: F' N! k: [3 g3 l% @提取码：z448


'''
" t# F; L& G& i% q, ?7 ]实验环境Python2.7.13,igraph包，cairo包，numpy包
) }9 n9 n+ O8 }'''
  ?7 i  m1 ]8 D! O# -*- coding:utf8 -*
5 ~" w1 g2 {: Ifrom igraph import *
import numpy as numpy
0 q" q: m  e2 q$ Z. X0 z  Kfrom  numpy import *
  v4 @5 d. u4 _  Zimport random
& {) m3 R* m+ t6 a) o+ o
def len_arr(infected_array,nodes_num):#获取感染数组长度
% J: }* r8 [. n9 `. q7 x    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

) F% z4 N0 ~/ H& v' q! j" bg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
5 x2 I8 q/ Z" c4 s2 B3 h5 u6 A. vsummary(g)
6 i% n* f0 }2 t1 g1 Znodes_num=g.vcount()#统计图中的结点个数
( [: O0 N% n) [. V7 x& e) Unet_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矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
5 L4 u% F5 Z* }( Lprint(nodes_state)
, _+ R) a* Q# q: Q; L2 }infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
% U. D5 j- ?- R6 h  E( a
: Q+ b) i. T( h1 I8 d* }: F' i' Pinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
4 e! N8 @9 _7 A, P  |' Zset_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
- l2 B, R* n  u! L  [nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
9 h& i0 z; B' T+ |, w9 \0 }nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
/ T: n* n2 |" ig.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

/ q( }- Z- L& q+ b) v* `stop=False#感染过程结束的标记
5 M" h. ]2 b# T8 a& k3 y1 h- m0 Ttemp_time=0#第几次感染
! F9 L0 M+ T' s7 X4 h- mtemp_len=0#本轮的感染源数量初始化
8 R" I% j! F* O9 t$ P
while not stop:
5 o% t. g- T0 R( N8 x8 d+ H    i=0#记录让每个感染源都传播一次
6 f1 l/ {( K$ D- U# c) D7 F+ H    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:
. l4 ?- R7 v2 F' B/ t+ {            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
2 v7 B$ ?) I" K) [4 g' x7 v& c            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
9 ?" T5 P- I& |& e6 B5 y            for j in range(nodes_num):#遍历节点
/ U! G0 n6 F+ U) S3 i5 S                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
* i1 J9 Z+ L( F                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
. d  h8 q( A3 y) z& Q            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:
                    nei_arr[t]=j+1
! C% v( w/ z' j9 a+ a                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
# @# H* w- v* J+ Y                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
' [- w* j1 }3 |  H: @            if len(ran_infe_arr)>0:#存在需要被感染的节点
. R5 d1 b6 L) U6 o                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
2 b# f4 y+ d: q) D+ k                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
2 U. U. o* F) J                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
8 d3 T# v; D; a" V  A                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
; i7 i' l3 q- Y9 t8 ]3 l7 A                    plot(g)#绘制
- h6 X" ?  @8 E0 b; M4 U& k" a                    t=t+1
4 o' [  \1 Y$ K/ e2 L& g            i=i+1
  P- W9 m2 s2 R; R# _    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
/ z3 a# e, L" K- X! J" A        stop=True


视频演示bilibili传送门
效果图
% W# w' X. p3 `: a! t' c3 r
8 B9 R; W9 W( E) b" X# |

4 `2 R) C# A# |% h0 s


1 `2 z" i; n$ }7 \: B& c  G

! [- O0 u+ M3 u& m/ A
/ S6 P; x) h- x
! G2 _: f* M" o: H% e% r
  x  }$ m+ v! Z; v7 O
. n- p2 r+ ^: p ————————————————
7 C" @, N7 W- R; j' A+ U! Y版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
% o" B7 K/ O& D: ]! C) L

尔雅

发表回复棒!
  E8 J7 h( E! j# x; P