Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
0 G  a% a9 B+ G6 r#SI疾病传播模型的原理
在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类
( h  l. r3 m+ _! [! v
易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
$ A( [- j: Q' l% R/ ZSI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
8 n* a/ `$ f# V即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
6 s8 w' w0 H* ~7 o在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
& G( d/ `& ]8 O9 j6 n8 d会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
( |' A& ~7 W9 K分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
显然有，
s(t)+i(t)恒等于1，此时之前的公式可以记做
ds/dt=-βsi
( f, Q: |, I* A  tdi/dt=βsi
! q7 z& |5 F3 k  _) ^即
di/dt=βi(1-i)
: z* q7 D: F: f. U6 U9 P上式也成为Logistic增长方程式(Logistic growth equation),
方程的解和图像如图

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


'''
0 e, i; {1 B4 R; F0 ^3 v实验环境Python2.7.13,igraph包，cairo包，numpy包
* j; {+ Q9 M" c3 F" {'''
! }5 n' T9 B* F+ Z$ J4 l6 P4 G# -*- coding:utf8 -*
from igraph import *
import numpy as numpy
from  numpy import *
7 i  {! v. w& H( J/ W% oimport random
6 A  b$ o: Y' |9 D. g# n$ c
  c2 B2 U  G* ]( q& @2 m  X& @def len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
  D- L: S; Q) X% J0 a, s" V    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
/ _9 s# y% H2 A- s( [    return len_value
$ V' v. B* G) T5 K
g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
, @$ k! ~8 Q) isummary(g)
9 [7 ^1 F: W' D5 K1 m) W* N' \nodes_num=g.vcount()#统计图中的结点个数
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
8 K) _/ |# @9 C2 i0 k4 s$ |g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
3 i6 f! k  `4 g: {a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
print(nodes_state)
3 X: ?* ^! V* a6 k9 k( B: minfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
6 r* P6 m* R  |1 X2 T/ }print(infected_array)
9 X3 D9 a% G+ _! W7 P( K
0 X6 G# \9 C4 {2 Winfe_rate=1#传播率(感染率) 1代表邻接点100%被感染
set_time=2#传播次数(感染次数) 2次
  F$ w3 Z3 r! P3 U, C6 ]4 J: o1 Vsource_seed=1#感染源位置
7 w% C' b) l/ w. _% l3 \% unodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
3 U6 k8 x' k7 r9 }! T3 Vnodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
9 ^: Z* T% t; D- i0 }4 {nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
& \: t+ p( I- L$ @$ _0 ~2 @g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
% b9 W" W3 `2 x) f/ rinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制

( P( ^- G$ y4 c3 Fstop=False#感染过程结束的标记
" |- ^; i! p( V" htemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化

& U) Q+ k- V7 C; ?9 R, H! t! @  C7 _- \while not stop:
    i=0#记录让每个感染源都传播一次
1 @$ A. e# y* k$ i8 @    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
3 }* N( r6 R$ f% u# [        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
, c7 F* }( {( \1 k2 K        while i<temp_len:
! S3 h) F' x$ m5 P5 N, T            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
# z4 Y* {9 {1 g% ~/ l/ I            nei_count=0#下一轮可以被感染到的节点数量
5 W2 I/ z5 q7 D6 s: O! A3 ^            #生成下一轮可能被感染的节点的集合nei_arr
  d4 b' W1 ~6 S1 r$ j            for j in range(nodes_num):#遍历节点
* J8 o/ i/ V* R( ^, e3 T. J/ U                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++
            nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
* E" I) x* h) u& O- |! a. L3 y            t=0
$ N* @/ b  [3 S/ ~, C            for j in range(nodes_num):
& N* T, A6 u5 I3 M2 X                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
" X" l2 p1 h1 W5 e6 o9 `) A1 [+ X                    nei_arr[t]=j+1
; L& Z6 C; N. B: _7 y                    t=t+1
            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
. W/ {2 C& C7 _: w$ D            if len(ran_infe_arr)>0:#存在需要被感染的节点
                t=0#让ran_infe_arr内每个感染源都被感染
7 B1 I% o6 @$ e, a" E- @2 V5 k                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
( g' O% h4 G; z4 ^# g, T                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
9 h! E& f+ L7 c% X* g7 d                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
, V; c- I4 p0 d# i7 p* z  j% ^                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
, L7 M" m: C! E4 y% E- ]                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
                    plot(g)#绘制
                    t=t+1
            i=i+1
* s( X' s% P% i) ~7 D& K    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
6 e+ O4 C0 z9 W5 H        stop=True
4 J7 L& ]1 i3 b- Y( g

视频演示bilibili传送门
效果图


+ t' _( Q& J5 h! W0 g
8 o8 ~# m- e  Z% b% [




, f/ A' l/ `2 b
* \* ~2 T3 y* A; w1 u9 n$ _  ?. ~

! I$ v# \4 @4 W/ N5 ^9 a
————————————————
3 N- Z# ?9 e8 o版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
" v( i  f5 F3 I, c* ^* J6 E+ O: O2 I
8 X$ r) Z; j$ f6 I
7 z2 g; s( g+ T) a5 D! }/ _

尔雅

发表回复棒!