격자 모델 시뮬레이션에서의 확산
시뮬레이션
이 포스트에서는 격자 공간에서 어떤 성분(Ingredient)의 확산 현상을 모방하려고 한다. 이는 그 동시에 SI 질병 확산 모델의 시뮬레이션이기도 하며, 공간이 제한되어 있다는 점에서 SIR 모델로도 볼 수 있다.
변수
- $t$: 현재 턴을 의미한다.
- $I(t) \in \mathbb{N}$: $t$ 턴에서 확산되고 있는 성분(Ingredient)의 양을 나타낸다.
- $S(t) \in \mathbb{N}$: $t$ 턴에서 성분이 없는 빈 공간(Space)의 양을 나타낸다.
파라미터
- $\beta \in [0,1]$: 확산률(Diffusion Rate)로써 성분이 확산되어가는 것을 나타낸다.
액션
(모든 칸에서, 매 턴마다)
- 계산: 각 칸을 기준으로 상하좌우에 성분이 몇개나 있는지 센다. $i$ 행 $j$ 열에서 계산된 수를 $n_{ij}$ 라 하자.
- 확산: 각 칸마다 $1-(1-\beta)^{n_{ij}}$ 의 확률로 성분을 가지게 된다.
코드 리뷰
Step 1. 초기화 및 격자공간 생성
julia> row\_size = 10
10
julia> column\_size = 10
10
julia> β = 0.5 # 확산률
0.5
julia> Random.seed!(0);
julia> stage\_lattice = rand(['S'], row\_size, column\_size)
10×10 Array{Char,2}:
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
julia> stage\_lattice[rand(1:row\_size), rand(1:column\_size)] = 'I'; stage\_lattice
10×10 Array{Char,2}:
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'I' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
격자 기반 시뮬레이션 튜토리얼과 비슷하게 격자 공간을 생성하고 랜덤한 위치에 성분을 넣는다. 위의 경우에에는 $6$행 $6$열에 생성된 것을 확인할 수 있다. 이를 히트맵으로 시각화하면 다음과 같다.
Step 2. 계산
julia> I\_lattice = (stage\_lattice .== 'I')
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
julia> count\_lattice =
vcat(I\_lattice[2:end, :], zeros(Int64, 1, column\_size)) +
vcat(zeros(Int64, 1, column\_size), I\_lattice[1:end-1, :]) +
hcat(I\_lattice[:, 2:end], zeros(Int64, row\_size, 1)) +
hcat(zeros(Int64, column\_size, 1), I\_lattice[:, 1:end-1])
10×10 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
주변에 성분이 얼마나 있는지 세어놓은 행렬 $(n_{ij})$ 을 만든다. 예제에선 벡터 연산을 이용해 구했지만 반드시 이렇게 계산할 필요는 없다.
Step 3. 확산
julia> probability\_lattice = 1 .- (1-β).^count\_lattice
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> hit\_lattice = probability\_lattice .> rand(row\_size, column\_size)
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
각 칸 별로 확산 받을 확률을 계산해서 행렬 $\left( 1-(1-\beta)^{n_{ij}} \right)$ 를 만든다. 이렇게 계산하는 이유는 주변에 성분이 많을수록 확산률이 높아지는 것을 반영하기 위해서다. 빈 공간에서 성분을 받을 확률은 $\beta$, 받지 않을 확률은 $(1-\beta)$ 다. 빈 공간 $i$행 $j$열이 주변 $n_{ij}$개의 성분에게서 확산을 받지 않을 확률은 $(1-\beta)^{n_{ij}}$ 이고, 결국 그들 중 하나에게서라도 확산을 받았을 확률이 바로 $\left( 1-(1-\beta)^{n_{ij}} \right)$ 이다.
julia> probability\_lattice = 1 .- (1-β).^count\_lattice
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> hit\_lattice = probability\_lattice .> rand(row\_size, column\_size)
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
이제 균등분포 $U(0,1)$에서 같은 크기의 행렬을 만들어 확산이 일어났는지 일어나지 않았는지 판정한다. 위 케이스에서는 상, 좌, 우로 성분이 퍼져나가는 것을 볼 수 있다. 실제로 히트맵으로 그려보면 다음과 같다.
이를 반복하면 본문과 상단과 같은 움짤을 얻을 수 있다.