SRC워크샵(ABM By Python)
===
###### tags: `Python`
#### <목차>
1. [Python에서 animation 구현하는 법](#animation)
2. [개체들간의 거리행렬 계산](#distance)
3. [배열](#array)
4. [SIR 모델에서 거리에 따른 감염확률 계산](#prob)
5. [SIR 모델 시뮬레이션](#simul)
6. [KDTree 최적화](#kdtree)

### 1. Python에서 animation 구현하는 법 <a id = 'animation'></a>

- 애니메이션 구현에 주로 사용되는 패키지
: from matploitlib.animation import ***FuncAnimation, ArtistAnimation***

- FuncAnimation
: 사용자 정의 함수로 그래프를 계속해서 업데이트 해나가는 구조
참고) [FuncAnimation 공식문서](https://matplotlib.org/stable/api/_as_gen/matplotlib.animation.FuncAnimation.html)
- ArtistAnimation
: 변화하는 그래프를 리스트를 인자로 주어 그려 나가는 구조
참고) [ArtistAnimation 공식문서](https://matplotlib.org/stable/api/_as_gen/matplotlib.animation.ArtistAnimation.html://)

def moving_agent(n):
    from matplotlib.animation import ArtistAnimation
    from IPython.display import HTML
    
    location = rand(2, n) # 위치 좌표 설정
    
    
    fig, ax = plt.subplots()
    ax.set(xlim = (-10, 10), ylim = (-10, 10)) # 범위 미리 고정
    
    
    loc_lst = [] # 변화하는 그래프를 담을 리스트 생성
    for frame in range(100):
        location += 0.2 * randn(2, n)
        loc = ax.scatter(location[0], location[1], color = 'blue', s = 10)
        loc_lst.append([loc]) # 변화하는 그래프를 리스트에 추가
        
    
    anim = ArtistAnimation(fig, loc_lst, interval = 5)
    return HTML(anim.to_jshtml())
        

### 2. 개체들간의 거리행렬 계산<a id = "distance"></a>

julia와 비슷하게 rand, randn 메서드를 사용할 수 있다.

from numpy.random import rand, randn

for문을 통해 거리행렬을 계산하는 과정, 코드 또한 julia와 일치하는 것을 확인할 수 있다

from numpy import zeros
from math import pi, sqrt

D = zeros([10, 10])
for i in range(n):
    for j in range(n):
        D[i, j] = sqrt(sum((location[:, i] - location[:, j]) ** 2))

D       

### 3. 배열 <a id = "array"></a>
이 부분 또한 julia에서의 코드와 거의 일치한다. julia에서의 push!가 python에서는 np.append()와 동일한 의미로 이를 대체하는 것을 제외하면 동일하다.

## 배열

state = np.array(['S' for _ in range(10)])
state = np.append(state, 'I')
state = np.append(state, 'I')


bit_S = state == 'S'
bit_I = state == 'I'
n_S = np.sum(bit_S)
n_I = np.sum(bit_I)

n = 12
location = rand(2, n)

D = zeros([n_S, n_I])
for i in range(n_S):
    for j in range(n_I):
        D[i, j] = sqrt(sum((location[:, i] - location[:, 10 + j]) ** 2))

D

for문을 돌리지 않고 distance메서드를 통해서 간단한 계산도 가능하다

from scipy.spatial import distance
pair = distance.cdist(location[:, bit_S].T, location[:, bit_I].T, 'euclidean')
contact = np.sum(pair < k, axis = 1)
print(pair)

for문을 이용해 계산하는 방법과 distance메서드에 의한 계산 결과를 비교해보면 동일함을 알 수 있다

- distance.cdist() : 두 배열간의 거리를 계산하기 위한 메서드
참고) [distance.cdist 관련 공식문서](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html)

### 4. 거리에 따른 감염확률 계산 <a id = "prob"></a>
contact 변수가 정상인이 감염자와 일정 거리안에 들어오는 횟수, 즉 접촉 횟수이다. 접촉하면 감염되는 확률을 $\beta$ = 0.5라고 할 때, 한 사람이 감염되는 확률을 $1 - (1-\beta) ^{contact}$로 표현가능하다.

b = 0.5 # beta

# 감염자가 되는 확률
1 - (1 - b) ** contact
infected = 1 - (1-b) ** contact < rand(n_S) # 난수를 발생시키고 이보다 작다면 감염자로 특정한다
infected

np.put(state, np.where(bit_S)[0][infected], 'I') # infected = True인 사람들의 상태를 'S' -> 'I'로 변환, np.where(bit_S) -> 튜플반환이므로 [0]으로 array만 뽑아옴
state

### 5. SIR 모델 Simulation <a id = "simul"></a>
거리가 변화함에 따라 정상인, 감염자들이 변화하는 SIR 모델을 정의한다
- S : Supseptible(정상인)
- I : Infected(감염자)
- R : Recovered(감염 후 회복 된사람)

def SIR(n):
    beta = 0.3
    sigma = 0.1
    mu = 0.1
    
    S_ = []
    I_ = []
    R_ = []
    t = 0
    
    ID = range(n)
    state = np.array(['S' for _ in range(n)])
    state[0] = 'I'
    location = rand(2, n)
    
    while sum(state == 'I') > 0:
        t += 1
        location += randn(2, n)
        bit_S = state == 'S'; n_S = np.sum(bit_S); ID_S = np.array(ID)[bit_S]; S_.append(n_S)
        bit_I = state == 'I'; n_I = np.sum(bit_I); ID_I = np.array(ID)[bit_I]; I_.append(n_I)
        bit_R = state == 'R'; n_R = np.sum(bit_R); ID_R = np.array(ID)[bit_R]; R_.append(n_R)
        
        D = distance.cdist(location[:, bit_S].T, location[:, bit_I].T, 'euclidean')
        contact = np.sum(D < k, axis = 1)
        infected = rand(n_S) < (1 - (1 - beta) ** contact)
        
        np.put(state, ID_S[infected], 'I')
        np.put(state, ID_I[rand(n_I) < mu], 'R')
        
    return S_, I_, R_

fig, ax = plt.subplots()

def animate(t):
    x = np.arange(t_end)
    y1 = S
    y2 = I
    y3 = R
        
    ax.plot(x[:t], y1[:t], color = 'blue', label = 'S')
    ax.plot(x[:t], y2[:t], color = 'red', label = 'I')
    ax.plot(x[:t], y3[:t], color = 'green', label = 'R')
    ax.set(xlim = (0, t_end), ylim = (0, 50000))
    
    
ani = FuncAnimation(fig, animate, np.arange(t_end), interval = 100)

HTML(ani.to_jshtml())

### 6. KDTree 최적화 <a id = "kdtree"></a>
KDTree는 데이터에서 가장 가까운 점을 탐색하는데 최적화된 트리 자료 구조로, 공간을 분할하여 자료를 저장한다. 이는 특정 점 (또는 그 주변의 점)을 탐색하는데 $O(log\ n)$의 시간복잡도가 소요되어, 일반적인 자료구조 $O(n)$보다 효율적이다.

1. ***KDTree.query()***
: 데이터 점들간의 거리 계산

## 최적화
import scipy
from scipy.spatial import KDTree

kdtree = KDTree(location[:, bit_I])
print(kdtree)
kdtree.query(location[:, bit_S].T, distance_upper_bound=k)[0] < k
# k보다 작은 거리에 True, 아니면 False 반환

2. ***cKDTree.queryballpoint***
: 데이터와 거리를 인자로 주면 거리보다 작은 데이터들을 list들의 array꼴로 반환

from scipy.spatial import cKDTree

tree = cKDTree(location[:, bit_I])
contact = tree.query_ball_point(location[:, bit_S].T, k, return_sorted=False)
for i, lst in enumerate(contact):
    if lst == [0]:
        contact[i] = 1
    else:
        contact[i] = 0
contact
    

SIR모델에 적용

def SIR2(n):
    beta = 0.3
    sigma = 0.1
    mu = 0.1
    
    S_ = []
    I_ = []
    R_ = []
    t = 0
    
    ID = range(n)
    state = np.array(['S' for _ in range(n)])
    state[1] = 'I'
    location = rand(2, n)
    
    while sum(state == 'I') > 0:
        t += 1
        location += randn(2, n)
        bit_S = state == 'S'; n_S = np.sum(bit_S); ID_S = np.array(ID)[bit_S]; S_.append(n_S)
        bit_I = state == 'I'; n_I = np.sum(bit_I); ID_I = np.array(ID)[bit_I]; I_.append(n_I)
        bit_R = state == 'R'; n_R = np.sum(bit_R); ID_R = np.array(ID)[bit_R]; R_.append(n_R)
        
        KDtree_I = KDTree(location[:, bit_I].T)
        contact = KDtree_I.query(location[:, bit_S].T, distance_upper_bound=k)[0] < k
        # contact = tree.query_ball_point(location[:, bit_S].T, k, return_sorted=False)
        # for i, lst in enumerate(contact):
        #    if lst == [0]:
        #       contact[i] = 1
        #   else:
        #      contact[i] = 0
        infected = rand(n_S) < (1 - (1 - beta) ** contact)
        
        np.put(state, ID_S[infected], 'I')
        np.put(state, ID_I[rand(n_I) < mu], 'R')        
    return S_, I_, R_

fig, ax = plt.subplots()

def animate(t):
    x = np.arange(t_end)
    y1 = S
    y2 = I
    y3 = R
        
    ax.plot(x[:t], y1[:t], color = 'blue', label = 'S')
    ax.plot(x[:t], y2[:t], color = 'red', label = 'I')
    ax.plot(x[:t], y3[:t], color = 'green', label = 'R')
    ax.set(xlim = (0, t_end), ylim = (0, 50000))
    
    
ani = FuncAnimation(fig, animate, np.arange(t_end), interval = 100)

HTML(ani.to_jshtml())

- 직접거리를 찾아주는 경우의 SIR모델 시뮬레이션 time
: <span style='background-color:#fff5b1'>6.2s</span>

- KDTree를 이용한 SIR모델 시뮬레이션 time
: <span style='background-color:#fff5b1'>1.8s</span>