Numpyを使ったプログラムをcupy化しようとしたら、
計算量多そうなscipyのcorrelate2dがあってどうなるかと思いましたが、chainerの畳み込み層を使ってみたらなんとかなった話です。

下の方にプログラム全体をまとめてあり、これを適当なpythonスクリプトとして保存して実行できるかと思います。

環境

Ubuntu 16.04
Python 3.5.2
numpy 1.14.3
chainer 4.1.0
CUDA toolkit 9.1
OpenCV 3.4.3

準備

六角形のグリッドでrock-paper-scissorsオートマトンを実装したこちら：
https://qiita.com/sage-git/items/183999c1006999959d68
これを普通の方眼グリッドでやってみます。

CPU版コード

まず普通のnumpyで組んだプログラムです。

import random
import numpy as np
from scipy import signal
import cv2

MAX_HP = 18

cell_color = [(255, 255, 255), (0, 0, 0), (36, 120, 1), (1, 84, 240)]

neighbor = np.array([
    [1, 1, 1],
    [1, 0, 1],
    [1, 1, 1]
], dtype=int)

class RockPaperScissorAuto(object):
    def __init__(self, W, H):
        self.W = W
        self.H = H
        self.F = np.zeros((self.H, self.W), dtype=int)
        self.HP = np.ones_like(self.F)
    
    def step(self):
        empty = self.F == 0
        rock = self.F == 1
        paper = self.F == 2
        scissor = self.F == 3
        nxp = signal.correlate2d(paper, neighbor, mode="same")*(rock + empty)
        nxs = signal.correlate2d(scissor, neighbor, mode="same")*(paper + empty)
        nxr = signal.correlate2d(rock, neighbor, mode="same")*(scissor + empty)
        self.HP -= (nxp + nxs + nxr)
        mut = self.HP < 1
        self.F[mut] += 1  
        self.F[self.F == 4] = 1
        self.HP[mut] = MAX_HP 
    
    def mutation(self, rate=0.0005):
        mut = (np.random.rand(self.H, self.W) < rate)
        self.F[mut == 1] += 1
        self.F[self.F == 4] = 1
        self.HP[mut == 1] = MAX_HP

    def make_image(self, cell_rad=10):
        width = self.W
        height = self.H
        img = np.zeros((height, width, 3), dtype=np.uint8)
        for i in range(4):
            img[self.F == i] = cell_color[i]
        return img

def main():
    F = RockPaperScissorAuto(480, 320)
    waiting = 10
    rate = 0.01
    while True:
        F.step()
        F.mutation(rate=rate)
        rate *= 0.99
        img = F.make_image(cell_rad=5)
        cv2.imshow("test", img)
        r = cv2.waitKey(waiting)
        if r == ord("+"):
            waiting = max(10, waiting//2)
        if r == ord("-"):
            waiting = min(1000, waiting*2)
        if r == ord("q"):
            break
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

六角形グリッドのときからの変更点

mutation

六角形グリッドの時と比べて、maskがなくなりました。

step

あるセルの周りのセルをカウントしてから、それに基づいて変異させます。

ふと、「HPが１以下になったものが変異する」とすると、mutationと同じルーチンで世代交代を表現できることに気づきました。なので後半3行はmutationとほとんど変わりません。この変更は前回の六角形グリッドのものでも有効です。

make_image

今回は丸を並べる必要がないので、セルをそのままピクセルに割り当てます。
imgのshapeが(H, W, 3)なのに対し、self.F == iのshapeは(H, W)です。これでもimg[self.F == i]というようにindexに入れることができて、余ったshapeの次元である3に対してcell_color[i]を代入しています。さすがnumpyのブロードキャスト。

cell_color

軍産複合体（参考：https://qiita.com/sage-git/items/183999c1006999959d68 のネタばれ）が導入した赤、青、黄の３色は目が痛くなったので、落ち着いた色に変えます。いろいろ探した結果、歌舞伎の定式幕になりました。

変異確率

mutationに渡すrateをステップごとに減衰させることで、初期に激しく揺動して渦を作り、後半は静かな環境で安定した渦を観察します。

これで実行してみた結果です。

典型的なBZ反応のシミュレーションができたかと思います。

cupy化

import

import numpy as np
import cupy as xp

openCVで画像表示をするので、numpyも使います。

initialize

np.zeros()など、numpyで作る配列をそのままxp.zeros()などのようにcupyに変えていくだけです。

correlateをconvolution2dで

Python：Scipyのcorrelate2D (convolve2d)でいろいろ数える
こちらで触れたように、correlateとconvolutionは場合によっては同じように使えます。
Chainerのconvolutionといえば畳み込み層。ということでconvolution_2dを使ってみます。
え？ニューラルネットのConv層は数学的には厳密には畳み込みではない？大丈夫です、数学的なことは問題になりませんでした。それにうまく行かなかったらこの記事を書いていません。

correlate2dで隣接マスの定義として使っていた行列をそのままConv層の$W$とします。
ただし、floatでないといけないので、とりあえず半精度浮動小数点にしました。
また、$W$の次元は$(c_I, c_O, h_K, w_K)$で、$c_I$と$c_O$は入出力チャンネル数です。今回はどちらも１なので、単純にreshapeで1を増やしておきます。

neighbor = xp.array([
    [1, 1, 1],
    [1, 0, 1],
    [1, 1, 1]
], dtype=xp.float16).reshape(1, 1, 3, 3)

まずrockなどセルの状態の数え方から。現在rockのセルが1、それ以外が0の行列を作ります。

    rock = xp.array(self.F == 1, dtype=xp.float16).reshape(1, 1, self.H, self.W)
    paper = xp.array(self.F == 1, dtype=xp.float16).reshape(1, 1, self.H, self.W)

これもfloatにします。また、convolution2dの入力は$(N, c, w, h)$ですが、バッチ数$N$とチャンネル数$c$は１としてreshapeします。

続いて、correlate2dと同値の関数を使って、あるセルの周りのpaperの数を数えます。

import chainer.functions as F
v = F.convolution_2d(paper, neighbor, pad=1)

なお、convolution_2dはchainer.Variableを返します。これは余計なデータを持っていますので、とっととint型のcupy配列になおしておきます。
ついでに、周りのpaperの数を気にするのはrockかemptyだけなので、それでマスクします。

    nxp = xp.array((v*(rock + empty)).array, dtype=int).reshape(self.H, self.W)

このnxpは、rockまたはemptyの周りのpaperの数で、この数がそのままHPを減らす量になります。
ここでrock, paper, scissorsに同様の処理をするのですが、このconvolution2dは並列にできそうです。ということはミニバッチ化できそうですね。

ミニバッチ化

もう一度数えるところから。

    rock = (self.F == 1).reshape(1, self.H, self.W)
    paper = (self.F == 2).reshape(1, self.H, self.W)
    scissor = (self.F == 3).reshape(1, self.H, self.W)

convolution2dの入力は$(N, c, w, h)$ですが、バッチ数$N$ 用の次元がreshapeから減りました。これを

    rps = xp.array(xp.stack((rock, paper, scissor)), dtype=xp.float16)

こうしてミニバッチ化します。stackによりミニバッチ用の次元が追加されます。ミニバッチ中のサンプル数は３です。これをconvolution_2dに渡し、cupyの配列に変換しつつミニバッチを分解します。

    cr, cp, cs = xp.vsplit(xp.array(F.convolution_2d(rps, neighbor, pad=1).array, dtype=int), [1, 2])

あとは必要な領域だけにmaskします。

    nxp = (cp*(rock + empty)).reshape(self.H, self.W)
    nxs = (cs*(paper + empty)).reshape(self.H, self.W)
    nxr = (cr*(scissor + empty)).reshape(self.H, self.W)

ここからの処理（HP減とセルの種類変化）はNumpy版と同じです。

image化

cupyの配列はそのままではcv2に渡せないので、GPUからCPUに渡してやります。

    img = xp.zeros((height, width, 3), dtype=xp.uint8)
    # 中略、imgに対して必要な描画を行う
    return np.array(cuda.to_cpu(img))

実行

プログラム全体

以上をまとめると、以下のようなプログラムになります。

cupy_grid_automaton.py

import random
import numpy as np
import cupy as xp
from chainer import cuda
import chainer.functions as F
import chainer.variable as V
import cv2

MAX_HP = 12

cell_color = [(255, 255, 255), (0, 0, 0), (36, 120, 1), (1, 84, 240)]

neighbor = xp.array([
    [1, 1, 1],
    [1, 0, 1],
    [1, 1, 1]
], dtype=xp.float16).reshape(1, 1, 3, 3)

class RockPaperScissorAuto(object):
    def __init__(self, W, H):
        self.W = W
        self.H = H
        self.F = xp.zeros((self.H, self.W), dtype=int)
        self.HP = xp.ones_like(self.F)
    
    def step(self):
        empty = (self.F == 0).reshape(1, self.H, self.W)
        rock = (self.F == 1).reshape(1, self.H, self.W)
        paper = (self.F == 2).reshape(1, self.H, self.W)
        scissor = (self.F == 3).reshape(1, self.H, self.W)
        rps = xp.array(xp.stack((rock, paper, scissor)), dtype=xp.float16)
        cr, cp, cs = xp.vsplit(xp.array(F.convolution_2d(rps, neighbor, pad=1).array, dtype=int), [1, 2])
        nxp = (cp*(rock + empty)).reshape(self.H, self.W)
        nxs = (cs*(paper + empty)).reshape(self.H, self.W)
        nxr = (cr*(scissor + empty)).reshape(self.H, self.W)

        self.HP -= (nxp + nxs + nxr)
        mut = self.HP < 1
        self.F[mut] += 1  
        self.F[self.F == 4] = 1
        self.HP[mut] = MAX_HP 
        
    def mutation(self, rate=0.0001):
        mut = (xp.random.rand(self.H, self.W) < rate)
        self.F[mut == 1] += 1
        self.F[self.F == 4] = 1
        self.HP[mut == 1] = MAX_HP

    def make_image(self, cell_rad=10):
        width = self.W
        height = self.H
        img = xp.zeros((height, width, 3), dtype=xp.uint8)
        for i in range(4):
            img[self.F == i] = cell_color[i]
        return np.array(cuda.to_cpu(img))
    

def main():
    F = RockPaperScissorAuto(640, 480)
    waiting = 10
    rate = 0.02
    while True:
        F.step()
        F.mutation(rate=rate)
        rate *= 0.99
        img = F.make_image(cell_rad=5)
        cv2.imshow("test", img)
        r = cv2.waitKey(waiting)
        if r == ord("+"):
            waiting = max(10, waiting//2)
        if r == ord("-"):
            waiting = min(1000, waiting*2)
        if r == ord("q"):
            break
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

パフォーマンス

系のサイズをいろいろ変えた時のGPUの使用率です。
GPUの使用率については、実行中にnvidia-smiコマンドを１０回程度打って、僕のインスピレーションによる平均を取ったものです。
GPUはGTX1060です。

右端のは確か2560x1960です。

これでnumpyのGPU化が必要であればcupyで、そして適宜chainerを使えばscipyのgpu化もできるのかもしれません。

Tensorflowでも同様にできる気がします。

その他

オートマトンの見た目ですが、一言、キモい。