Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
# K4 i% v* \/ T$ u#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
. a' F1 ^. c% N0 N; M; F
* V7 y9 N0 O4 F& m易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
0 C5 ^$ q# v5 m" g( s. A$ Q易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
" G# w* y& V% I4 ?" \" w( J& BSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
3 n6 y# C  O- H3 k+ h$ d2 p, D即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
* S8 r; U0 ~# F! c7 T& Z( Hds/dt = -β*S(t)I(t)/N
/ Y* @6 ^' u3 d# r3 c同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
; C! j# k7 k& u4 K/ ]$ e分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
7 C" {+ V7 F8 B  g/ m- \* A! p显然有，
) i7 w4 k. z4 H5 b( _) O: I4 Ms(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
8 x" c! x5 o8 Y9 P/ ?; y; \: W8 jdi/dt=βsi
即
7 [) E. j! A# u5 K# r; p% Qdi/dt=βi(1-i)
4 j( q& U3 n. |上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

4 S/ |4 S3 w( o" T* D6 L代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448


'''
/ A5 d2 G; {0 j' s实验环境Python2.7.13,igraph包，cairo包，numpy包
* P5 x; F1 X& g; d  e2 C, u( C'''
# -*- coding:utf8 -*
from igraph import *
" z0 j* y/ a- b. S* m- G/ L+ ?import numpy as numpy
+ r+ {+ Z# l9 I" n1 B* y- Afrom  numpy import *
, n- p7 [' N+ \$ X+ ]  J& B2 Aimport random
. {' {/ C6 s& y& l. p
def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
; M" P) _  Q4 F# b2 y# ]! w    return len_value
8 h/ [& c% A1 b/ W
4 y0 L  _$ c5 z- r5 }/ R) yg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
, W) z7 y/ p0 p2 hsummary(g)
nodes_num=g.vcount()#统计图中的结点个数
/ W* X) d2 u. Fnet_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矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
5 c$ G, C. W/ ^                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
' d' J% l$ Q& F8 Y- jinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)

/ F: q3 t* h$ iinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
- k. r6 [" H9 q6 `# _3 Sset_time=2#传播次数(感染次数) 2次
+ Z0 Y  J6 M; \source_seed=1#感染源位置
' v1 W0 N9 ?& ^" K+ Znodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
, \: p! F. ~. D( t$ A: }2 v; ]nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
- x+ K7 w& T8 @6 d# {8 C3 H+ i2 }* Bnodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
) I8 O/ X2 R# e" x( Tinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
. \# X) |  d/ t8 y$ J) y0 W( Qplot(g)#绘制

stop=False#感染过程结束的标记
temp_time=0#第几次感染
4 e2 H4 Z) }9 Q/ i% b1 ~9 ntemp_len=0#本轮的感染源数量初始化
. v3 U& r' P+ R- U2 c& j5 q
7 f- L/ G/ _. A& o) b) L: }4 Cwhile not stop:
1 P  T! }; Z2 u8 p7 Z  u/ W1 [    i=0#记录让每个感染源都传播一次
    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)#获取本轮的感染源数量
  T/ B& S- |" G/ {1 i3 L        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
4 e; n2 k8 S; R                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
# ]3 e* L: c) y: v8 l/ p, |- ]7 R            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
% b6 ^% k& u- Z7 _( J1 G6 |            t=0
            for j in range(nodes_num):
. K; X% P+ n  V! D* p0 z) T$ u                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
1 V* I2 G0 u" N6 E' y                    nei_arr[t]=j+1
                    t=t+1
# e3 e5 S" }8 I2 n            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
) }, Y7 X* w. t/ n3 i, W5 @                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
( J6 ?! }6 P9 ]" r                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
$ c  R2 m" k6 X, N$ X( \. {7 E                    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"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
( e. f. R( b7 }4 C- n& h& `6 c            i=i+1
    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
% g, w( i6 r7 n: ~: i

% c* W$ Y  P. B5 |( F2 j; D视频演示bilibili传送门
- e8 T  `. H: R效果图


8 `7 _9 M% l1 F" R6 F2 N
) f2 k0 l6 K" J, Z3 D




* g2 p, K) E" b' J* N

4 Y' ~$ }+ \8 \
' O* s; f4 F* M/ u+ X
————————————————
2 S1 o- _" b# g4 t版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
8 ]- ?& a$ R; L# v1 z. B原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

% t0 s! H8 L' y# s1 }. d6 M

尔雅

发表回复棒!
' j: Y0 x6 Z1 f9 q% h# F! z