Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
9 `: D8 k+ T& L# H#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
* f# W% Y8 v: G& T9 t+ ?移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
# f9 @; Q+ [* D, k$ NSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
- W, {. K: m* p# S6 J  e8 [在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
  B# x# h2 ?( C9 I% v$ Q1 u, c6 cS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
) L/ S6 O+ N6 b2 {# x  D! E/ d感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
& b% \$ Q& E( Vds/dt = -β*S(t)I(t)/N
# s  `# |; _* v" e5 k同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
3 c. M; |" j7 a: P: C分别将时刻t处于易感状态和感染状态的个体所占比例记为，
# ^4 y, [- a5 Gs(t)=S(t)/N
i(t)=I(t)/N
显然有，
6 X0 [: O, D) R( W" zs(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
di/dt=βsi
. A6 C7 G6 K# v4 h6 M- F/ s即
5 a& g; R' K: Z7 ?' R( k& Pdi/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
, ?: e& o; \) O" h9 M' s5 r方程的解和图像如图

代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
) d" H* i* p/ r. E+ f提取码：z448
/ A- w0 I$ L8 A1 h
& \0 ~6 W" H' A& D# c
'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
7 ^$ ~6 M! ~9 ]( ]9 z'''
# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
7 t$ ~$ q, B1 kfrom  numpy import *
import random
! q2 n( ~  n0 y0 S& T: {1 c
def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
: X: e' @8 ^' q* R+ o
0 j- `6 L" p, Q- H- |g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
+ i9 k$ }$ K9 G; O- gsummary(g)
& s) ?" S! M* j# ]# l! mnodes_num=g.vcount()#统计图中的结点个数
net_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表示，第三列是节点感染的时间
, Z1 Z* ?, V+ j$ |# R7 E( Kprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
! r) {6 f. N7 E0 h0 \print(infected_array)

infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
9 @' T  t8 K% e$ q$ _- I  o9 M, hset_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
, B9 z( H! J" ~nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
) ?; m: \' I* U- t6 [$ k$ ]nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
1 L: |; |! i2 e* `2 X* `infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

5 n! W- |6 z: Z& N" Q# Vstop=False#感染过程结束的标记
% i) }* P; L: Gtemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
- u! j6 Y/ _  |8 @+ t# s1 N6 z    i=0#记录让每个感染源都传播一次
) G1 O* y. T2 Y) c) R' {  {    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
9 K" x6 Y) b. w, d" r        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
) z  L. b  [. U6 l# p: |7 a        while i<temp_len:
/ K) n% I4 ?9 u+ c& ?4 ~* B            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
1 M( J8 l% J$ _' E            nei_count=0#下一轮可以被感染到的节点数量
9 `% o) G9 d) @1 V8 Q' ]            #生成下一轮可能被感染的节点的集合nei_arr
, B; b3 a, c; F            for j in range(nodes_num):#遍历节点
5 J3 `; A! l0 X; {! C0 _0 ^                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
/ I; h! \. u9 u' U, i) N) U                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
& X. v; I# w8 K4 L1 Z: ^- W            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
& T) x: p( }/ k7 h" F1 W$ e" H            t=0
, w7 T: D  R' M# h. \/ e            for j in range(nodes_num):
2 i4 a4 A/ K4 g7 G5 T                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
7 M: n1 |0 Q' H' P5 S7 ?* f# E                    nei_arr[t]=j+1
1 j& x9 ^' J4 U1 o0 m; r                    t=t+1
! H- P2 \% f; z! M6 I& [" I            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
% M: Z! l% W8 N            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
+ j1 U3 F$ O) G8 s) O% D' c0 ?4 B                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
) z* k0 X! R# Z! O, I0 a                    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"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
- v2 |2 f( O: D            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
) A% g8 [5 K& R  O9 g        stop=True
5 p8 m/ |6 D5 @1 v0 W
4 ?% c# }; w- W
视频演示bilibili传送门
# Z7 r8 I+ u9 B& B$ j* A9 {效果图

5 W! F6 A  b& N; e0 R
% S. l8 B7 n; G: i

+ s6 {) V- k0 ]0 X
# V" S' O6 t) x; d



8 \8 \2 o6 P( y) O8 i! Z" e


————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
; K- ~' p% \6 _) Q( d* O8 c原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
+ I& c- W$ c. B

5 C2 @1 ^. M: p* O; f* U0 _8 s+ A

尔雅

发表回复棒!