- 在线时间
- 1630 小时
- 最后登录
- 2024-1-29
- 注册时间
- 2017-5-16
- 听众数
- 82
- 收听数
- 1
- 能力
- 120 分
- 体力
- 563412 点
- 威望
- 12 点
- 阅读权限
- 255
- 积分
- 174246
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 5313
- 主题
- 5273
- 精华
- 3
- 分享
- 0
- 好友
- 163
TA的每日心情 | 开心 2021-8-11 17:59 |
|---|
签到天数: 17 天 [LV.4]偶尔看看III 网络挑战赛参赛者 网络挑战赛参赛者 - 自我介绍
- 本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。
 群组: 2018美赛大象算法课程 群组: 2018美赛护航培训课程 群组: 2019年 数学中国站长建 群组: 2019年数据分析师课程 群组: 2018年大象老师国赛优 |
Python实现简单的SI传播模型
+ s6 { A5 t* Y#SI疾病传播模型的原理
- ]# l4 X: }" O; G在经典的传染病模型中,种群(Population)内N个个体的状态可分为如下几类
3 C; B2 o9 }+ M& F3 l4 h# U9 d& G8 ~3 \- l
- M9 ?/ s! g* j* h! u易感状态(Susceptible)。一个个体在感染前是处于易感状态的,即该个体有可能被邻居个体感染。* \! S; e/ Y1 C9 n
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。,即该个体还会以一定概率感染其邻居个体。
" N# b9 \/ \ W- x4 a$ |' B6 j移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态,当一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑改革提。5 Q F, v Q( T
SI传播模型是最简单的疾病传播模型,模型中的所有个体都只可能处于两个状态中的一个
% O9 j0 `5 x1 b; _3 M1 j: \即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。( k" p3 W1 }4 D' S3 N
在给定时刻t,令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目,显然有
P9 |9 P/ ^8 J: ^6 iS(t)+I(t)恒等于N,这里,N是个体总数。随着时间t的增长,易感个体与感染个体的接触# f7 T u5 M, r
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β,疾病仅在
" J- B4 R+ _6 G3 J' m6 s, O感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t,易感个体的比例为S(t)/N,感染个体的数量为I(t),一次,易感个体的数量将以如下变化率减少
; h# R/ t" }, f: M- W9 Gds/dt = -β*S(t)I(t)/N9 r. ^1 }* F' w! T3 }) {
同时,感染个体的数量会以与易感个体相反的变化率增加,) H5 O. t: S M; O* M+ r
ds/dt = βS(t)*I(t)/N
2 ^7 ~( r2 E# ]; ~7 o) h# q分别将时刻t处于易感状态和感染状态的个体所占比例记为,
$ k; @ x6 L: x7 m# @* E+ v2 ys(t)=S(t)/N
1 k A) x% q! ji(t)=I(t)/N
5 }+ O' B* @9 _显然有,
( t" s% x2 F& Z3 S$ \( H) q8 Is(t)+i(t)恒等于1,此时之前的公式可以记做
% k& }8 X4 S! t; a& U% a A% h9 ]; Wds/dt=-βsi9 k" W# o+ A9 H) r; t% q/ ?, N
di/dt=βsi% d9 |% e- i! R6 w8 R
即
( `' P8 ~( y9 t# H c* F' Adi/dt=βi(1-i)
+ I" t- b: ]8 k$ q5 h8 @# w- Y上式也成为Logistic增长方程式(Logistic growth equation),
! `" J2 `! k' w' u5 w3 G方程的解和图像如图
- C( g8 w/ [; |$ r- m7 B, y4 j
5 A7 V% l8 j& {1 A代码和相关文件以及环境链接:链接:https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
+ c! E" s A' h. |5 \提取码:z448
" T, }6 k ~/ o: f/ v/ i
: v4 h$ u" G( e; {5 v; d) B% P" G$ ?5 p5 n( ^% t( V4 j
'''7 g9 `% L& d0 g! S
实验环境Python2.7.13,igraph包,cairo包,numpy包& k$ k& c, f+ w! E0 b4 G
'''
* V' G `- n/ N0 M4 A# u4 j$ ]# T# -*- coding:utf8 -*
/ N4 S+ J$ T; h) D. K" `; h n& rfrom igraph import *' _& N" R3 C* l3 O+ g% c
import numpy as numpy8 V s; U$ V* x8 \: i$ Z
from numpy import *
& Z9 ]# L( L0 Pimport random
O/ K/ ~9 b3 {9 i! O0 \8 K E, J
def len_arr(infected_array,nodes_num):#获取感染数组长度3 C4 E" }& a& h. b+ V( m) ] H
len_value=0#初始化长度- U+ f6 c5 ~/ f3 G( ^$ Q% R
len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)4 D& j9 J+ d( M8 J' b7 P$ U
return len_value4 K+ d, l1 r: j& d# _7 W
6 L% |( H# A* F- Dg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
3 `7 G) v6 A, W: s( Zsummary(g)
3 e P. O- |+ g! i' k* x- {nodes_num=g.vcount()#统计图中的结点个数
4 r# Q, z- R) N* i5 snet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
: g% T1 @, } g4 N% a! b/ Qg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
% K1 W5 Y$ ^, Ba=[arange(nodes_num)+1]*3#声明一个N行3列的数组a; C f6 L7 ?; o
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建,用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重!!!)
0 r _: C5 ^" l #第一列是节点编号,第二列是节点状态,感染状态用-2表示,第三列是节点感染的时间
# c3 P, {% H* R" H8 h7 Hprint(nodes_state)
. y5 c3 z+ R: d/ [; n" Rinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染 34代表网络节点数
1 U, i5 B; J/ m4 F! cprint(infected_array)
2 c4 p) i2 F; t4 c9 x
& q( w2 ~2 Q( n( ]8 _$ x4 Uinfe_rate=1#传播率(感染率) 1代表邻接点100%被感染; ]: ~. k; \! M$ A3 p3 V2 |
set_time=2#传播次数(感染次数) 2次
" o! D" c$ o; }7 c' ssource_seed=1#感染源位置
; g5 i P# V; Z) {: X+ Hnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
& {" I/ l7 O0 V2 n" inodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
# S* P/ L0 m# u, @# t/ o" Unodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
9 x7 U, D4 y$ X' S# U# u% }. s3 @g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
5 i6 a- a6 C9 |% }& }infected_array[0]=source_seed#将感染源的位置存入被感染节点列表! d- Z* ]# ?) T t2 I9 R
plot(g)#绘制
9 Q+ N6 ]: x5 o3 e" A2 t3 a w5 o2 g6 Z% @# }
stop=False#感染过程结束的标记3 L; B6 D0 L" J* V" w
temp_time=0#第几次感染$ e& K3 |; r* q) @* S% N7 w0 _
temp_len=0#本轮的感染源数量初始化
5 V: E7 F/ p. Q" F% M |1 I& h4 C5 t( R2 m
while not stop:" f3 C, w1 `) o- `
i=0#记录让每个感染源都传播一次
+ U# z$ i/ n# f! N: W5 B: }% |, N if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行" k0 Z' L/ F; h1 Y- @* a3 y+ B
temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
% Z M5 S/ J8 a0 s/ o while i<temp_len:2 W8 P5 G/ Q) ~ Y5 Z) L' o
temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间7 n0 |, t8 w( v
nei_count=0#下一轮可以被感染到的节点数量! E6 U7 r1 Z. p0 j+ q8 f/ D/ D
#生成下一轮可能被感染的节点的集合nei_arr
- v# \3 ^1 J( W/ P( Y+ ` for j in range(nodes_num):#遍历节点
; Q4 z: d! J4 L* O' C if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染+ {- l3 E" z2 q$ G$ V! `' X
nei_count=nei_count+1#下一轮可以被感染到的节点数量++
! A# I# B$ o, f4 N P# j nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染+ }8 q$ c8 l* }
t=0
+ \' g8 t, l" z' M$ \$ K6 I' t for j in range(nodes_num):
4 c# F8 z2 I q/ @- P- Y if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:9 a$ |- r( `1 q0 ^
nei_arr[t]=j+1
0 |2 R5 Q7 D# ]5 X4 F t=t+1
& ~( O9 x$ q" b ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
- n! K2 b1 M& ~ #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
- j( U4 h$ X/ l& i3 q( V% } if len(ran_infe_arr)>0:#存在需要被感染的节点
5 K+ B0 u Z" q5 M& [# n t=0#让ran_infe_arr内每个感染源都被感染
9 ~( U9 \7 e; `3 | while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
) r5 x* s. P, E9 Z nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态 m- R4 J( `. t' \, E8 w1 z
nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间. H/ S) l% h. f% A
infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中8 q: b. E9 }& O `- l2 V
g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
' x5 u/ [" M$ x. l1 ]) E$ [! ^ a plot(g)#绘制
* ^' A7 r( w% C% o, b t=t+1
$ n/ E% I3 ?7 N& Z | i=i+1
+ R m+ {# O4 T$ ^ if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染3 e m2 G/ s3 I3 V) {/ d
stop=True
3 ~! ], D6 B' s$ d! i' a* y3 j/ \/ z
! y& m' _: V) ]% }% X3 y2 J3 u$ a
视频演示bilibili传送门; l5 E H! Z1 F. {0 j
效果图2 h( A2 s* Y4 Z
+ _! S) M6 ?+ B( m& p( X& V- L
, Y9 R4 B# D9 M7 {. ]
$ n. ?& T: A3 v. Z, u/ `
2 n- R F. b9 l$ B( @
5 o' O6 \( P+ K1 [5 C0 P
5 H: w0 T* z' w" P7 q% N
" y; c: f$ i E3 O- r7 J
4 b- [- M6 B9 M+ [% b; S
% Q4 ]/ l. g7 v) w. y
" q% G6 K q1 L
; N0 Q! c5 _1 Z3 b( T: U# B5 R
4 P0 y& j$ K$ S& |* q' h
————————————————
( g* ?1 D2 l+ C" t1 K9 S版权声明:本文为CSDN博主「eck_燃」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。5 y/ g) `' S3 ^% N4 f) U7 w
原文链接:https://blog.csdn.net/wdays83892469/article/details/80878862. W" [6 z0 @) H D7 B
6 `6 E/ V, f- i! h- B
( U6 _! F' N& F6 J |
zan
|