Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
/ j9 h5 D# T( N1 K#SI疾病传播模型的原理
6 u9 v; s/ |. J: g1 b1 b' S1 D7 ?在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
; {% q8 S, P  E) s6 Y  o% ^
易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
; u+ [/ g' B4 ?移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
( |: L( n: y. I) @8 ]% D即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
( }1 X, B/ p  P, A在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
* e/ L! E( ?* P7 }# BS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
, ]6 ?4 @/ }% c# _  G: s' `1 N6 v感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
  m% F- q. v) ^; Q3 K6 bds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
0 @, i. z6 I/ Z1 g$ u. `ds/dt = βS(t)*I(t)/N
: n1 v1 y6 C1 V& e4 S分别将时刻t处于易感状态和感染状态的个体所占比例记为，
9 ^) d! W" i. |& k; n# u9 Ns(t)=S(t)/N
" [0 @) T1 o8 s6 U0 o! ?7 Di(t)=I(t)/N
显然有，
, V/ `- P* Y, s0 D: ~4 h: s  Vs(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
9 _6 b! ^! y" o. i3 edi/dt=βsi
9 o& N; N8 C8 \8 e2 @即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
9 W, D2 D  U4 E: x) e: g# g( X方程的解和图像如图

  z5 `( x# }- _% |9 ?! D0 N代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448
( B" \! F- ]: d, k- A' @

'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
- b# i; t9 o, R& j: s9 ^'''
# -*- coding:utf8 -*
0 L$ Z0 F! X0 ~" t( |; c) |+ jfrom igraph import *
import numpy as numpy
from  numpy import *
import random
1 `1 s8 c% O# l2 v7 {* l/ s9 A
& D- I3 C& a0 T8 `def len_arr(infected_array,nodes_num):#获取感染数组长度
' J7 ^0 `# c) @5 A4 p4 d+ I) L. |    len_value=0#初始化长度
+ b* ]/ P/ V/ L& M( n# L8 v7 L- E    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
  G& C! N* K" T5 w8 V
. |6 N: e5 g6 v( l% `4 Gg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
1 }4 ]# x- \7 M: N* k9 Enodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
4 o# z# Q+ Z' ea=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
; D$ m# m4 |1 b5 H" D; K9 l4 ~6 G                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
8 ]0 G0 r9 {9 u. K' Uinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)
+ J3 T* u. E0 e. g6 D  j
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
- r. y0 p5 I3 m. X: c8 B9 S( L8 hset_time=2#传播次数(感染次数) 2次
/ @! q9 Y5 }- g& Q# m2 lsource_seed=1#感染源位置
4 v) l0 s7 C% s$ K! Mnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
5 k0 k' v: a' y/ s% R6 _nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
/ q* \: r  q+ T2 t# Kg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
# B' K& s) n. Z5 ]plot(g)#绘制

; F) y% e$ ~! tstop=False#感染过程结束的标记
temp_time=0#第几次感染
0 q9 C: @" @+ J( ]7 ^' w7 k6 n5 ftemp_len=0#本轮的感染源数量初始化

while not stop:
    i=0#记录让每个感染源都传播一次
/ M8 c9 V# o! k7 R/ ?    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:
- I5 a( o& y, N" o. g            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
5 C8 l  D: }! \: |: j( K! u            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
6 O' s+ g2 p+ c1 }5 U+ ~            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#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
& q) ^- y- _% E+ s# p  k' q            t=0
" x) R8 R4 A8 K& W. Q            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
$ ~5 c! I! t2 J7 _2 u. Z$ U) H                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
/ H) B" h/ [* C* l                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
) ?, v/ @/ |& i* H            if len(ran_infe_arr)>0:#存在需要被感染的节点
( R4 r' Q1 z5 c% w) J                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]]#将此次感染节点放入总的感染节点数组中
6 w# n1 U  }) G                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
3 H- ^" W) a' N0 z& K            i=i+1
+ C% n& M; e5 e6 ^- @    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
# _/ p4 y; G7 I# A) N! I& Y8 `        stop=True


视频演示bilibili传送门
8 ~  A# N# M: |/ x2 N/ F效果图
; r0 w7 v4 [$ z/ p1 m$ K* A. K2 l9 H
6 K) R) r- z% L
- Q( d' |' p% E* h; k: T* A




  E. u, z2 @( U2 h


; |& S0 r: X7 k: |+ _7 h: Y
, o* m" \% n" M8 n( L
————————————————
% O/ n" ^! G- G/ Z版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
" J7 L! U* O0 z0 W& c原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
0 R7 y/ b& h) i  n$ s6 n
) j0 |  L& Q2 J" R
5 f7 t% a/ l" x# y

尔雅

发表回复棒!