Python + Numba CUDA でライフゲーム

Python + Numba CUDA で GPU プログラミングしてみるシリーズ

- その０： [Python + Numba CUDA で GPU プログラミング入門以前](https://qiita.com/hoto17296/items/332da1d57173534c3aa7)
- その１： [Python + Numba CUDA で画像処理](https://qiita.com/hoto17296/items/d7f3bd3df0a51bf78bf9)
- その２： Python + Numba CUDA でライフゲーム **[← イマココ]**

前回までで一応 GPU で計算はできるようになったけど、あの程度のことは CPU でも余裕でできてしまうので、今回は **GPU のありがたみが感じられる何か**を作ろうと思った。

GPU を使うアプリの題材としてライフゲームが面白いと思ったので、今回はこれを実装する。

## ライフゲーム
[ライフゲーム \- Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%A4%E3%83%95%E3%82%B2%E3%83%BC%E3%83%A0) より引用

> **ライフゲーム (Conway's Game of Life)** は1970年にイギリスの数学者ジョン・ホートン・コンウェイ (John Horton Conway) が考案した**生命の誕生、進化、淘汰などのプロセスを簡易的なモデルで再現したシミュレーションゲーム**である。単純なルールでその模様の変化を楽しめるため、パズルの要素を持っている。
> 
> 生物集団においては、過疎でも過密でも個体の生存に適さないという個体群生態学的な側面を背景に持つ。セル・オートマトンのもっともよく知られた例でもある。

> 我々が暮らす空間、さらには時間が連続的なものであるか、それとも非連続的なものであるのか、という問いはギリシア時代から思索の対象となってきた。セル・オートマトンはその問いに答えるものではないが、**空間、時間が不連続であった場合、どのような世界が形成されるのかを示してくれる。**

ライフゲーム、エンジニアなら一度はハマる (と思っている)。

## 実装
GUI の実装は Python 標準モジュールの Tkinter を使った。

ソースコードを全部載せると長いので詳細は以下のリポジトリを参照。
[hoto17296/python-gameoflife](https://github.com/hoto17296/python-gameoflife)

CUDA の処理の部分だけ抜粋する。

``` python:gpu.py
import numpy as np
from numba import cuda
from base import GameOfLifeBase


@cuda.jit
def update_cell(before, after):
    """対象のセルの次の状態を計算するカーネル関数"""
    x, y = cuda.grid(2)
    # 周囲のセルの生存数をカウントする
    count = 0
    for i in (-1, 0, 1):
        if x + i < 0 or x + i >= before.shape[0]:
            continue
        for j in (-1, 0, 1):
            if i == 0 and j == 0:
                continue
            if y + j < 0 or y + j >= before.shape[1]:
                continue
            if before[x + i, y + j]:
                count += 1
    # 対象セルの状態を周囲のセルの生存数から、対象セルの次の状態を決定する
    after[x, y] = count in (2, 3) if before[x, y] else count == 3


class GameOfLifeGPU(GameOfLifeBase):

    def __init__(self, *args, **kwargs):
        # 親クラスの init の時点で self.state の初期状態が設定される
        super().__init__(*args, **kwargs)
        # 現在の状態と次の状態を入れるための配列を GPU メモリに転送する
        self.device_state = (
            cuda.to_device(self.state),
            cuda.to_device(np.zeros(self.size, dtype=np.bool)),
        )

    def update(self):
        """状態 (self.state) を更新するメソッド"""
        # カーネル関数を実行する
        update_cell[self.size, 1](*self.device_state)
        # 最新の状態をデバイスからホストにコピーする
        self.state = self.device_state[1].copy_to_host()
        # 次に update 時にカーネル関数の引数に渡す順番を入れ替える
        self.device_state = tuple(reversed(self.device_state))


if __name__ == '__main__':
    app = GameOfLifeGPU()
    app.start()
```

前回の画像処理と違う点としては、 `cuda.to_device` を使うことで ホスト → デバイス の状態の転送は最初の一回だけで済むようにした、つもりなのだけど・・・。

実際に試してみると、毎回ホストからデバイスに状態を送った場合と処理時間に差がなかったので、GPU メモリの扱い方を勘違いしているか、もしくはデータ転送以外の処理がボトルネックになっていると思われる。

ひとまず動くものは出来上がった。

<blockquote class="twitter-tweet" data-lang="en"><p lang="ja" dir="ltr">Python + Numba CUDA でライフゲームが実装できたぞ ＼＼\\٩( &#39;ω&#39; )و//／／ <a href="https://t.co/C99ydrUW6i">pic.twitter.com/C99ydrUW6i</a></p>&mdash; ほと (@hoto17296) <a href="https://twitter.com/hoto17296/status/992075801439883264?ref_src=twsrc%5Etfw">May 3, 2018</a></blockquote>

## CPU 版の実装
比較するために CPU だけで計算するバージョンも実装した。

``` python:cpu.py
import numpy as np
from scipy import signal
from base import GameOfLifeBase


kernel = np.array([
    [1, 1, 1],
    [1, 0, 1],
    [1, 1, 1],
])


class GameOfLifeCPU(GameOfLifeBase):

    def update(self):
        """状態 (self.state) を更新するメソッド"""
        # 各セルの周囲のセルの生存数をカウントする
        count = signal.convolve2d(self.state.astype(int), kernel, mode='same')
        # 各セルの次の状態を計算する
        self.state = self.state & (count == 2) | (count == 3)


if __name__ == '__main__':
    app = GameOfLifeCPU()
    app.start()
```

**シンプルすぎわろた。**

状態 (bool) を数値 (int) に変換した行列 `self.state.astype(int)` は各要素が 0 または 1 であるため、行列

```math
\begin{bmatrix}
1 & 1 & 1 \\
1 & 0 & 1 \\
1 & 1 & 1
\end{bmatrix}
```

を畳み込むだけで「周囲のセルの生存数のカウント」が全セル分まとめて計算できる、というのがポイント。

こんな簡単に書けてしまうとやっぱ CPU って楽だな〜と思ってしまう。
(NumPy と SciPy がスゴいという話ではある)

## CPU との比較
空間サイズを 2000x2000 (400万) セルに設定して実行したときに「1秒あたりに何回状態を更新できたか」を計測した。

| 実行環境 | 計測結果 |
|---|---|
| GPU あり | 約 13.0 回/秒 |
| CPU のみ | 約 4.9 回/秒 |

GPU の方が速い！すごい！

けど・・・、あれ？
思ったより差が小さいな・・・。

CPU 版はシングルスレッドなので CPU 1コアしか使ってなくて、それでもこの程度の差しか出ないのは不思議。

おそらく、以下のような理由だと思われる。

- NumPy と SciPy がスゴくて CPU の性能を最大限に引き出せている
- GPU 版の実装にチューニングの余地がある

## まとめ
Python + Numba CUDA でライフゲームを実装したら、CPU と比べて 2.6 倍くらい速かった。
まだまだ改善の余地がありそうだけど今回はここまで。