セル・オートマトンによる感染シミュレーション

セル・オートマトンは格子状のセルと単純な規則による離散的計算モデルで、あるセルの時間tにおける状態、および近傍セルの時間tにおける状態により、あるセルの時間t+1(次世代)における状態が決まります。ルールを確定的ではなく確率的に変えた確率的セルオートマンを用いることで感染症の広がりをシミュレーションできます。すなわち「隣のセルが感染者の場合の感染確率をxとする」といったルールを設けることでt→t+1(現世代→次世代)での感染の広まりを計算することができます。

具体的な計算手順は
- 現世代の状態を記録する配列変数を用意し初期状態を記録する。
- 初期状態の全セルと隣接セルの状態を確認し、それを元に次世代の状態を計算する。
- 計算した結果を次世代の状態を記録する配列変数に記録する。
- 全セルの次世代状態を記録し終わったら、現世代の状態を記録する配列変数にコピーする。
となります。
セルの数は縦横100セルの1000セルとし、初期状態では中央の1つのセルを感染者(セルの値=1)、その他のセルを未感染者(セルの値=0)としています。全セルについて隣り合うセルに感染者がいれば20%の感染率で感染します。上下左右それぞれについて感染者をの存在をチェックし、感染者がいれば都度感染判定を行っているため、隣り合う感染者が多いほど感染確率は高くなります。また、端の取り扱いについては配列変数を102x102と一回り大きくつくり、周囲のセルの値は常に0とし処理します。
計算結果は200世代分の図をアニメーションgifにして感染が広がる様子を視覚化しています。
なお、実際の新型コロナウィルスの感染状況を反映してはいません。

関連記事

SIRモデルによる感染シミュレーション SIRモデルによる感染シミュレーション(一度低下した感染率が再上昇する場合)

環境

• windows10 home
• Anaconda 3/ jupyter notebook 5.6.0
• Python 3.7.0

コード

import numpy as np
import random
from PIL import Image, ImageDraw

cells = 100  # 縦(横)のセル数
side = 2  # 1つのセルの一辺ピクセル数

cell = np.zeros((cells+2, cells+2))  # 各セルの状態を記録するndarray
cell[cells//2+1, cells//2+1] = 1  # 中央のセルを1(感染)に設定

trans = np.copy(cell)  # 各セルの次のステップの状態を記録するndarrayにコピー

rate = 0.2  # 感染率
time = 200  # 時間

# 　各セルの状態を描画する関数
def image():
    im = Image.new('RGB', (cells*side, cells*side), (0, 0, 0))
    draw = ImageDraw.Draw(im)
    for y in range(1, cells+1):
        for x in range(1, cells+1):
            if cell[y, x] == 0:  # セル値0(未感染者)は青
                color = (0, 0, 255)
            if cell[y, x] == 1:  # セル値1(感染者)は赤
                color = (255, 0, 0)
            draw.rectangle(((x-1)*side, (y-1)*side, x*side, y*side), fill=color)        
    image_list.append(im)

image_list = []  # 各時間ごとのセルの状態画像を保存するリスト　
image()  # 各セルの状態を描画する関数を呼び出し

for t in range(time):
    for y in range(1, cells+1):
        for x in range(1, cells+1):
            if cell[y+1, x] == 1:   # 注目セルの下のセルをチェック    
                probability = random.uniform(0, 1)
                if probability < rate:
                    trans[y, x] = 1
            if cell[y-1, x] == 1:   # 注目セルの上のセルをチェック     
                probability = random.uniform(0, 1)
                if probability < rate:
                    trans[y, x] = 1
            if cell[y, x+1] == 1:   # 注目セルの右のセルをチェック     
                probability = random.uniform(0, 1)
                if probability < rate:
                    trans[y, x] = 1
            if cell[y, x-1] == 1:   # 注目セルの左のセルをチェック     
                probability = random.uniform(0, 1)
                if probability < rate:
                    trans[y, x] = 1
    cell = np.copy(trans)  # transをcellにコピー
    image()  # 各セルの状態を描画する関数を呼び出し

# アニメーションgifファイルに保存
image_list[0].save('cellular_automaton.gif', save_all=True, append_images=image_list[1:],
                   optimize=False, duation=200, loop=0)

実行結果

未感染者(セルの値0)を青、感染者(セルの値1)を赤で表示しています。感染の空間的な広がりを視覚的にとらえることができます。

続いて、感染者が回復するモデルを計算します。感染者(セルの値1)が回復者(セルの値2)に変化する割合として回復率(0.05)を導入します。回復者は再び感染することはありません。

コード

import numpy as np
import random
from PIL import Image, ImageDraw

cells = 100  # 縦(横)のセル数
side = 2  # 1つのセルの一辺ピクセル数

cell = np.zeros((cells+2, cells+2))  # 各セルの状態を記録するndarray
cell[cells//2+1, cells//2+1] = 1  # 中央のセルを1(感染)に設定

trans = np.copy(cell)  # 各セルの次のステップの状態を記録するndarrayにコピー

rate = 0.2  # 感染率
recover = 0.05  # 回復率
time = 200  # 時間

# 　各セルの状態を描画する関数
def image():
    im = Image.new('RGB', (cells*side, cells*side), (0, 0, 0))
    draw = ImageDraw.Draw(im)
    for y in range(1, cells+1):
        for x in range(1, cells+1):
            if cell[y, x] == 0:  # セル値0(未感染者)は青
                color = (0, 0, 255)
            if cell[y, x] == 1:  # セル値1(感染者)は赤
                color = (255, 0, 0)
            if cell[y, x] == 2:  # セル値1(回復者)は緑
                color = (0, 255, 0)

            draw.rectangle(((x-1)*side, (y-1)*side, x*side, y*side), fill=color)        
    image_list.append(im)

image_list = []  # 各時間ごとのセルの状態画像を保存するリスト　
image()  # 各セルの状態を描画する関数を呼び出し

for t in range(time):
    for y in range(1, cells+1):
        for x in range(1, cells+1):
            if cell[y, x] != 2:  # 回復者は除外する
                if cell[y+1, x] == 1:   # 注目セルの下のセルをチェック    
                    probability = random.uniform(0, 1)
                    if probability < rate:
                        trans[y, x] = 1
                if cell[y-1, x] == 1:   # 注目セルの上のセルをチェック     
                    probability = random.uniform(0, 1)
                    if probability < rate:
                        trans[y, x] = 1
                if cell[y, x+1] == 1:   # 注目セルの右のセルをチェック     
                    probability = random.uniform(0, 1)
                    if probability < rate:
                        trans[y, x] = 1
                if cell[y, x-1] == 1:   # 注目セルの左のセルをチェック     
                    probability = random.uniform(0, 1)
                    if probability < rate:
                        trans[y, x] = 1
            if cell[y, x] == 1:  # 回復判定
                    probability = random.uniform(0, 1)
                    if probability < recover:
                        trans[y, x] = 2
    cell = np.copy(trans)  # transをcellにコピー
    image()  # 各セルの状態を描画する関数を呼び出し

# アニメーションgifファイルに保存
image_list[0].save('cellular_automaton(r).gif', save_all=True, append_images=image_list[1:],
                   optimize=False, duation=200, loop=0)  

実行結果

未感染者(セルの値0)を青、感染者(セルの値1)を赤、回復者(セルの値2)を緑で表示しています。

以下のサイトを参考にさせていただきました

セル・オートマトン『ウィキペディア（Wikipedia）』
セル・オートマトンによるパンデミックモデルの構築（その1）
Pythonの文法メモ >> 【標準ライブラリ】randomによる乱数生成
Pythonの文法メモ >> 【Pillow】ImageDrawモジュールによる直線、四角、円、多角形描画

ブログランキングに参加しています

にほんブログ村