Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
+ b5 Q- V. H3 m. R, _& a( |) s在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
: t- Q" M9 ~. E; a8 h
& g* j' y! t9 J易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
; A9 h$ y" {6 |+ G4 u: U, I( aSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
- i  N3 g2 w* G  g; z) W  y即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
2 ~$ [1 c- K* s+ r% m# hS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
+ X4 P! M7 u( I# G同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
- D; e0 O6 q; l2 O! r) x' g分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
. x9 c2 U& x6 R6 qi(t)=I(t)/N
# O1 `  J) Y, P6 n+ a. f5 |显然有，
$ m0 |6 f; L; y! Ts(t)+i(t)恒等于1，此时之前的公式可以记做
7 N8 o6 Y) S2 }) dds/dt=-βsi
di/dt=βsi
) X0 P" g8 X! J" f2 d即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
3 @; Z/ Y. n6 J2 j+ _9 u# H& m5 D方程的解和图像如图
% o' I3 e+ g4 A7 w/ V6 L
& J+ A3 g; V4 |, H) B! }代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448
0 e- N0 L5 |6 Q2 A' n* x
  u4 K( C' y$ j. @  [4 R
'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
7 Q3 m1 H% f2 i. Q) `  Z3 x, {'''
# -*- coding:utf8 -*
7 m. o5 h4 G" g; n/ @from igraph import *
  a& q' g" r5 q/ v' o0 kimport numpy as numpy
2 u) d- [) E3 q+ Rfrom  numpy import *
import random
" U7 w  \8 c* B: ^/ p4 ?
6 _* N9 U9 [) _, tdef len_arr(infected_array,nodes_num):#获取感染数组长度
: a8 `6 Z& t# @0 C* n3 z6 _    len_value=0#初始化长度
" G1 o! S$ h# m! A    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

. P$ P! Z1 F5 Mg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
( b- |9 t/ Q/ d1 w* h4 Tsummary(g)
  ^" F/ l0 X2 L7 _/ S8 e$ Pnodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
1 u" w7 k/ T$ fg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
& w0 U" y; ]" z2 |a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
/ ]; P8 R4 n' X3 o                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
' S% h  Y7 C/ T$ f% B4 H6 V$ v1 }, m6 Uinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
, f4 f% _- D+ M0 D1 Yprint(infected_array)

1 H: x0 d4 q$ c  p/ \) v# P/ dinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
' z) I+ ]9 M8 xnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
+ S$ A& ~4 r' o7 _nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
0 {9 {0 h: S: ^( K) B  Gg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
: D) E& B+ X5 I! ?/ iinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

0 ~1 U! n; U: w" d! H& \% |stop=False#感染过程结束的标记
# D1 ~  l4 G" Y/ Jtemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
0 x$ ?' |1 Z( M+ _, y9 Q    i=0#记录让每个感染源都传播一次
$ [% o9 S0 Y, ?. D! I    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
$ |' v4 n$ A; @) ^' G2 ~        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
        while i<temp_len:
+ Q) p% V7 Y8 Y3 W7 Z            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
1 ^' w( j5 h" j/ T            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
  D! J. s: v& R# ?            for j in range(nodes_num):#遍历节点
: X$ _6 Z. w& _% ~/ W7 t/ a                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
* w! a# c2 n7 E( B                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            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
                    t=t+1
5 Y% O+ Q% i2 F            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):#对刚才生成的会被感染的数组内的节点进行感染
                    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"#将此次感染的节点集的所有节点颜色置为粉色
- ?* E& ^5 y! J; T6 C& i                    plot(g)#绘制
                    t=t+1
            i=i+1
$ C2 s# a- c# o, t/ _* w    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
5 B5 o8 `2 Y" d$ y        stop=True
* q3 m7 T" t! [2 W: \7 \, Z, L
+ i) v  \- f( H# U
视频演示bilibili传送门
" q. K. c6 t& G' r  J效果图
5 Z. Z2 s& G0 L3 K; [: f







; H4 _8 m2 H- I# B+ s* |
* t8 e2 j# Q; _: {9 b


————————————————
* Q/ n$ t# D: u) |版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
7 N" _  U: I9 ]原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

- Y: r) ^7 K4 _  V( `
( N9 d" h/ K3 [. m/ O! U; C& t7 O

尔雅

发表回复棒!
1 a2 A" b) \" W; E- U5 f4 T