Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
/ d2 W1 J8 X* C4 H1 e- u( N#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
4 P2 c7 K. }8 p( Q2 B
/ j, D9 b5 o  y/ G. Q' v4 A4 M' t易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
. t$ E$ q8 `$ J# X% p移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
7 f0 g" \! v2 g9 Y7 Z即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
' O5 f4 d: t  P$ u) P9 c( z会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
4 E2 P/ Z! O7 P  y8 Z; _ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
8 E2 y" y" _. `ds/dt = βS(t)*I(t)/N
: U- h/ ^' o. r. {) M7 \% u  w分别将时刻t处于易感状态和感染状态的个体所占比例记为，
+ g5 G# `  J* Os(t)=S(t)/N
, L3 n; N/ X/ J! P# }% _$ ^i(t)=I(t)/N
显然有，
5 S  y5 h: b9 c/ f  n5 M3 \5 rs(t)+i(t)恒等于1，此时之前的公式可以记做
- P+ V: m$ @8 H; e0 F" v3 Hds/dt=-βsi
di/dt=βsi
% q( P  }+ P/ }# J- Q( t% y即
di/dt=βi(1-i)
& b$ m3 Y" C* e3 l. ~; q1 f上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
8 f8 u) F0 A# S4 F! f提取码：z448

  J" d  r+ F, o6 X
% Q# A2 F& T) n5 z'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
; h& [- j, ]" m$ u! i. J3 I'''
  `! R+ F6 L  Y; K% D/ T, q# -*- coding:utf8 -*
from igraph import *
; i' b. m: X! D! S# B& x% U6 V* zimport numpy as numpy
! ?1 a1 o4 _: S: M) ?5 Efrom  numpy import *
import random
6 i+ f) t0 W  B# a
def len_arr(infected_array,nodes_num):#获取感染数组长度
1 F* ~: P) d  f7 x3 {/ m6 j    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
% U$ |  O  q4 t    return len_value

' x9 V: x( ^9 @2 a: S$ T" Ng=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
$ F3 r" x1 n5 Z% h/ N& @" \# H5 Ynodes_num=g.vcount()#统计图中的结点个数
2 k4 Z2 J- `- v3 r: V. L- Pnet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
- P9 x1 K& ]" \' F& J$ L" }1 qg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
0 r# N2 q" ]. o3 u7 \" y0 Wnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
9 P5 T' b+ k' ?2 p: f* H( Z' Oprint(infected_array)

) b" r) U; K9 U8 J" rinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
7 L% ]/ m" `5 ^/ i9 ^* Wnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
" S  H$ ^8 e$ j; x# C* }nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
+ m2 i' I9 u# T9 s) m; Rg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
5 P9 Z' o4 T. O% |0 R* pinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
$ U4 O7 }; `" e# Hplot(g)#绘制
% W( n, W* z# c0 p0 F  u, r
1 d' a* ?# l- `2 |( b7 a- xstop=False#感染过程结束的标记
temp_time=0#第几次感染
; }; q1 ?3 I  Z. J( Z) }temp_len=0#本轮的感染源数量初始化

' K. G7 |. P; m7 G1 D0 Y% S  R8 j0 qwhile not stop:
7 V# E, Z5 z$ m2 u    i=0#记录让每个感染源都传播一次
( P$ w$ |5 F# |% |1 l' G& [. n    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)#获取本轮的感染源数量
' s  m/ \6 T* J- U        while i<temp_len:
' W3 w) ~# v' s# y8 z, l            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
) `4 I, }$ I, F* z" ?            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
& _. O( p" e2 Z0 V            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:
4 E; Y, S7 }9 \8 i5 ?                    nei_arr[t]=j+1
) m* y( S7 k" b: r                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
( H  {; f3 g- R% x                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
$ G( e; D9 k% S7 A  ]            if len(ran_infe_arr)>0:#存在需要被感染的节点
0 t$ @3 P( \' c# Y                t=0#让ran_infe_arr内每个感染源都被感染
9 y' h& L  Q% A( u/ s                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
                    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]]#将此次感染节点放入总的感染节点数组中
                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
( X; f+ [4 {) [5 ?/ X: d                    plot(g)#绘制
                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
0 R* P) K/ m+ S( o4 H
: M' r0 t. [# p# c( H
视频演示bilibili传送门
效果图
) E1 ^* |. e4 _
( ~, b9 j( @/ _3 x
1 O! }: l( ?+ u2 y


% c( c$ I. Q5 c9 a
, ^4 O4 T$ y5 _$ o



$ f8 U2 ~. ~+ v; M( U6 d4 _
# s1 C, s: X- `, _; P
————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

2 e7 r$ r) L9 p5 h2 R

尔雅

发表回复棒!
8 N* u3 j) O* w, c. g; V% s