Python实现简单的SI传播模型

杨利霞

Python实现简单的SI传播模型
2 M5 {% t7 M4 z# S; `7 D#SI疾病传播模型的原理
8 ?6 X# U3 t% Z在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类

8 ?1 F% F" V( S. I. {0 k1 k$ M$ ^# T易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。
0 @1 w+ v) v; D- X2 ~, x* USI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
9 M% B. o& Y7 l/ N% d0 x; A3 m" L在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有
S(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在
% k  _2 Y1 r* g1 T感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少
ds/dt = -β*S(t)I(t)/N
, k  J! v4 {( S* ]+ x; l同时，感染个体的数量会以与易感个体相反的变化率增加，
ds/dt = βS(t)*I(t)/N
分别将时刻t处于易感状态和感染状态的个体所占比例记为，
s(t)=S(t)/N
i(t)=I(t)/N
/ w) Z: I+ S. d显然有，
( E% [/ I7 Q& y* q( D( j& b5 |s(t)+i(t)恒等于1，此时之前的公式可以记做
; M; v% x/ a5 I# L) t0 ?ds/dt=-βsi
  M- w+ y0 F4 V9 ]" x/ W2 Xdi/dt=βsi
即
di/dt=βi(1-i)
上式也成为Logistic增长方程式(Logistic growth equation),
+ ~+ @/ a3 \5 k" ?6 k; d方程的解和图像如图
4 d9 D% N" c- U/ |( L! c9 t/ K
8 a# u( T6 b2 j7 r/ A代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
提取码：z448

4 \3 J  F: \0 ~- D- u
' l( T( H# L8 Q'''
$ F/ Z3 |4 ]& B5 i: |! k  K8 D实验环境Python2.7.13,igraph包，cairo包，numpy包
'''
8 n' @( z& K* _! P. E5 [- N# -*- coding:utf8 -*
from igraph import *
( L" O6 o* j. a$ yimport numpy as numpy
from  numpy import *
import random

6 H, t6 D- m- k) g" l# q- C8 Sdef len_arr(infected_array,nodes_num):#获取感染数组长度
    len_value=0#初始化长度
    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    return len_value
1 E, K) h& p; K3 v' c+ ^
# j/ ^9 ~9 K* ]g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
summary(g)
nodes_num=g.vcount()#统计图中的结点个数
, u$ I- B; {& @6 c% Onet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
! ]0 N& S( c6 N' h5 g% @3 `; Rnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)
/ E* F3 x" A- b$ E                                                    #第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间
5 v9 w2 |) L) ^1 b1 A1 a% pprint(nodes_state)
8 Y* N3 H+ g2 zinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数
print(infected_array)

infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
: [/ E) Z6 b2 W; ^* ]2 T" ^  hset_time=2#传播次数(感染次数) 2次
' c- q8 _' B) xsource_seed=1#感染源位置
9 ~3 x( r, |; w& T$ }% Anodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
/ Z$ x) B$ Q  _3 N& S! g5 ninfected_array[0]=source_seed#将感染源的位置存入被感染节点列表
plot(g)#绘制
" ?3 ~8 W0 ^1 Q6 U5 y0 b: _  J6 m
stop=False#感染过程结束的标记
$ A" E; K; u) `: ^' S3 f2 Ktemp_time=0#第几次感染
temp_len=0#本轮的感染源数量初始化
. J# f: o& B! [) v, h8 O
% w" \1 t4 C1 ~while not stop:
0 `+ u: v4 o: z" |; {# U( v    i=0#记录让每个感染源都传播一次
) m: L' b5 P. I8 n/ S' j6 S7 A    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
5 V* m/ ], z- I5 U- l- P        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
  G: A4 H6 o' X3 P0 @        while i<temp_len:
            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
( G. l( ?6 a! E3 V1 D, S            nei_count=0#下一轮可以被感染到的节点数量
' M5 L8 d5 l  O7 E6 M7 {  [            #生成下一轮可能被感染的节点的集合nei_arr
            for j in range(nodes_num):#遍历节点
- d  \0 \1 b4 j& }  c+ s3 r                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
            for j in range(nodes_num):
                if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
4 F. u" Z" K( p- `( h. P                    nei_arr[t]=j+1
) B* y3 k. `7 A1 B; O4 y6 F                    t=t+1
7 i" Z* b- y9 T' E            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
3 a. g  U1 g8 K3 p- @8 X                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
" f2 ]- R& ~8 t8 Y* Y( L            if len(ran_infe_arr)>0:#存在需要被感染的节点
9 i8 E) \7 X  ], b. M                t=0#让ran_infe_arr内每个感染源都被感染
                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
( V2 c$ u" Z3 ?" ]                    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]]#将此次感染节点放入总的感染节点数组中
                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
' S  h7 k/ @" @                    plot(g)#绘制
' y3 J( i0 K  Y$ P                    t=t+1
            i=i+1
8 _7 D0 x: U" P" ~% a    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
/ I* H6 f, F# U4 G# u+ C        stop=True
9 E2 w3 Q) J- e, Q$ M: f  [& C" S

" a7 n/ m( ~2 H2 e, H( o$ X* s/ W视频演示bilibili传送门
7 [3 k$ _% k/ `4 ^$ V7 E$ b& c效果图
. i* F9 V9 V: p% v4 l

# o8 K, _+ E8 U( B/ Z


9 q( n( d. z- X5 ^, _, A" ~" k1 N
( x: e8 o' A8 U2 Y

% ^0 e* O) R2 }7 ~* c* H
6 P: x5 S3 O( n4 F% q) H7 C: a( X  P


- O" S, N+ A- U3 {! i ————————————————
1 s1 \5 z0 @$ f% D! u版权声明：本文为CSDN博主「eck_燃」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862
! s9 e7 Z+ }: p. z3 L  m& O5 G

尔雅

发表回复棒!