Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
* \8 Y( Z0 S$ c" ^( T% z#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

, _7 V% E1 O9 O' Z$ @( N2 ]易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
+ f7 w2 G9 w! ]+ X9 E" f易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
6 l0 y8 E- |! P! oS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
- D1 w, Z, _6 S  ^) p9 p感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
6 m9 }. Y: p' q: z8 ?ds/dt = βS(t)*I(t)/N
/ U- l! R0 p2 h/ `+ D分别将时刻t处于易感状态和感染状态的个体所占比例记为，
% y, [/ g! n2 B3 z, Vs(t)=S(t)/N
i(t)=I(t)/N
显然有，
, Z9 H4 Z" F8 d5 c. O+ Vs(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
4 _  F% j$ x3 y+ ?di/dt=βsi
即
di/dt=βi(1-i)
/ v5 l" U. V" W2 {$ w; e' R上式也成为Logistic增长方程式(Logistic growth equation),
( X9 Z# H! I, S. o; J方程的解和图像如图

4 n3 F2 {8 D9 `- W# }代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448
3 [# {1 r5 k$ b. J7 j& R: _
5 L, p4 {3 A! E' X
'''
. o7 z. ~0 c, D& ?7 v2 i) U; o, f) ^实验环境Python2.7.13,igraph包，cairo包，numpy包
8 M) B/ M* l: r; k- L/ r'''
' f0 x2 V: q- Z# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
" {4 B0 D, C) Pfrom  numpy import *
+ ?) A1 |$ L" \; Gimport random
: h2 P& ]# E: O9 x9 ?; N2 E8 ]# r
def len_arr(infected_array,nodes_num):#获取感染数组长度
: L1 Y9 }* _: S. y    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value

  w( m$ F. x2 |5 R5 g, mg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
9 |( {1 E, z/ T, F$ Sg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
. \8 g: W; c+ X: B7 W/ S$ i8 w; Da=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
) Z; T4 E9 z0 ^( r! V7 knodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
# i$ m( I- v! l2 O; |: N- a! @$ g7 L                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
, n# y/ i2 d8 W$ P( e
% K, ^' R; o/ Sinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
  l! c3 O: l& I2 e% b/ e7 zsource_seed=1#感染源位置
& L( s; N0 c" z* T9 d  ?! Nnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
7 w) S: J6 `  b6 G$ vg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
  u9 ?9 r0 i3 h, L  N1 a5 Hinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
6 T& r  g# A0 A4 r( o0 lplot(g)#绘制
& b) K  B. ]( g
stop=False#感染过程结束的标记
temp_time=0#第几次感染
! @7 v* P1 _8 I+ M! m' [# h. L8 }temp_len=0#本轮的感染源数量初始化
- m. z& B9 ]) ^8 O- m
while not stop:
6 S9 x4 u4 G6 a8 F1 m8 J) L0 w    i=0#记录让每个感染源都传播一次
    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
# H! R- E% T* @; n* A        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
4 e4 r) s, ?; r% l8 T            nei_count=0#下一轮可以被感染到的节点数量
+ s, d7 `! m) y+ b            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
, `5 c9 e3 i" B" \& I                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
            t=0
$ K! \3 D4 o! q3 e1 s            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
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
$ `( T7 j0 W! P( w0 S1 d) P8 Z                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
            if len(ran_infe_arr)>0:#存在需要被感染的节点
$ |1 D9 u* l! W0 S3 M                t=0#让ran_infe_arr内每个感染源都被感染
. C* g% d! b9 x: V" S) A: \; p( B                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#记录感染时间
4 Y$ A5 e# k8 T- s1 v  N) s                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
# v! ^/ ~  Y: J                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
* i' \8 c( d, W( A  ]5 h                    t=t+1
            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
+ n4 I" U& j! ~, d7 d9 ]$ b        stop=True
/ R" M6 ^! x+ q$ [
- p; ]( d2 y! q; `  b$ q. Q3 @
视频演示bilibili传送门
效果图

! E# Q# R. n  O2 H" y8 S
  i8 g! A6 ^: X7 d5 y
9 d: u+ h, P9 ~: B5 l9 m
$ x6 m) T; H' n9 _! r

0 _; Z3 V  Q/ y) i

9 `  c4 R' a- R* l' @5 _# w4 U2 ^1 W
+ Y4 D, p5 U. ~1 l" b, t


————————————————
版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
6 r: V1 U# T7 M2 S原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

尔雅

发表回复棒!
; X6 O: E. W% V( }& p" ~% h. E