Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
#SI疾病传播模型的原理
! A) g+ l5 L! p$ p在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

/ l2 X* @) R) `4 L/ `易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
0 a6 Y& N% [+ V: @会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
* p5 l2 i$ C  D& ?ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
$ l% P; K2 Y1 K% A" Cds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
( t4 Z3 C, {8 ^2 w$ Ss(t)=S(t)/N
: @( [% g- K, u# pi(t)=I(t)/N
+ ^. |4 _6 v2 r显然有，
+ ?2 k! _) @# d% E" is(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
+ B3 G& W; @* r; m5 U$ ndi/dt=βsi
% W) c( L6 M; |8 `, U7 Y即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
9 N* `4 Q& Y- e% `0 j3 C5 G方程的解和图像如图

代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448

& n2 x- h" V$ r+ b; x: P
'''
实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
# -*- coding:utf8 -*
! q  p, A* s" C2 L- t2 a. _  a2 T. Ifrom igraph import *
import numpy as numpy
from  numpy import *
" z% Z' v1 V: D' B0 e* Q7 simport random
; l5 F, j) w  `: c
def len_arr(infected_array,nodes_num):#获取感染数组长度
% F& R$ y/ ?, O) f% O5 n. J1 a0 z    len_value=0#初始化长度
6 C  `, @8 w3 Q  Y+ f# z9 H    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
& a" R1 U: ]# n+ |. ~  B& K6 o' o    return len_value
& z/ I- r3 S/ q- e6 u1 T9 W
' d8 [- v9 H  A) U( ng=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
# t- F! i1 L# s% u) mnodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
! {" F5 C% d. ^0 ^3 T* g+ `2 Fg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
$ g) e* Z4 @5 {nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
! A9 ]  Z: `+ _/ W9 a* E  X                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
' e8 b! e+ x# H- O4 e7 Bprint(nodes_state)
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
9 m, x$ L& r' v, Xprint(infected_array)
  L  j, v( \$ y0 p% j
% K% S0 }4 `- `8 |infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
8 t5 e- D/ {  y5 |set_time=2#传播次数(感染次数) 2次
source_seed=1#感染源位置
) |2 a  x  y; d9 d, C9 hnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
9 p- `  k: G. a& h9 U$ w& Z2 Enodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
3 N$ l- g  I4 z% h1 z" ig.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
1 [6 H, B& @* P8 h1 ~1 P0 cplot(g)#绘制

stop=False#感染过程结束的标记
- C1 b6 D( ?5 b( Z: Atemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

while not stop:
    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)#获取本轮的感染源数量
        while i<temp_len:
' D) M% ]3 R$ Q7 i+ a% N            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
            nei_count=0#下一轮可以被感染到的节点数量
5 s, ~( A8 z8 F+ U7 i* j0 i            #生成下一轮可能被感染的节点的集合nei_arr
. b6 f, P" y# Q/ O: c            for j in range(nodes_num):#遍历节点
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
" c$ R1 C$ V  h6 U# U% N                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
7 i6 Q) P  V* ^/ A' w. ]            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
- X6 ]* u% S; s( A# h* }            t=0
            for j in range(nodes_num):
( b: Q/ i7 F6 Q  I                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
                    nei_arr[t]=j+1
2 o4 J/ x3 n" `* v                    t=t+1
3 C% a( m) c8 h6 Q5 y/ n9 W            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
4 V/ c2 J& V, d* v                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
* d/ s6 x9 p( ~3 o4 o( T3 z            if len(ran_infe_arr)>0:#存在需要被感染的节点
6 k. r: {/ [9 w* E. P! m4 K" o                t=0#让ran_infe_arr内每个感染源都被感染
7 n9 w7 r/ R$ ]. [" A; u9 E                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#记录感染时间
3 p; c; {. V# g8 y$ W4 W3 U& B                    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)#绘制
% {$ L& {6 T& d1 o% ?8 s: F                    t=t+1
- W/ Q7 W( e0 A  C" W+ b            i=i+1
8 a( J9 ?: a) ^2 F( ]/ Y    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
        stop=True
9 Z% }1 O4 F) V
/ `. @" E4 [6 \' R0 ^4 g- a' I6 X
3 F' U0 H3 e2 h5 }$ o! d0 i视频演示bilibili传送门
效果图

  L  j- P$ ^  b1 o* f- E
! `" ~$ I; h) J5 a' W- v5 n( B! K

2 b  D+ ]& B" @# W- h- ~
( v" z/ A6 t# z
- r- l5 k/ F6 h0 ]) L  ~9 K$ p


9 w9 w& ?; J% A3 _5 c


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

/ t. Q/ j0 Y5 W% Q' p

尔雅

发表回复棒!
1 q9 D6 d& P% D7 A) U; d- v, v