Python：Scipyのcorrelate2D (convolve2d)でいろいろ数える

Numpyには畳み込みの計算をするconvolve関数があります。ですがこれは１次元のみにしか対応していません。
一方でScipyにはcorrelate2D, convolve2Dが提供されています。

この定義をうまく使えば、前回のpythonでライフゲームで、あるマスの周囲９マスの１の数を

 signal.correlate2d(F, mask, mode="same", boundary="wrap") 

と一行で書けるというのが分かったので、今回はライフゲーム以外のケースでも試してみることにします。

HEXマップ

普通の四角形の方眼について、1行毎にマス目の半分だけずらすと、六角形格子と同じ並びになります。
六角形格子に対してどのように座標を振るかは諸説ありますが、今回は斜めにずれていくタイプとします。

    height = 8
    width = 10
    grid = np.array(np.random.rand(height, width) + 0.2, dtype=int)

このように定義した8x10サイズの格子を、

def print_hex(grid):
    H, W = grid.shape
    for i, v in enumerate(grid):
        print(" "*i + ("{} "*W).format(*v))

と定義したサブルーチンで出力してやると、

1 0 1 0 0 0 0 0 0 0 
 0 1 0 0 0 0 0 0 0 0 
  0 1 1 0 0 0 1 0 0 1 
   0 0 0 0 1 0 0 1 0 1 
    0 1 0 0 0 0 0 0 0 0 
     1 0 0 0 0 0 1 0 0 0 
      1 1 0 0 1 1 0 0 0 0 
       0 1 0 0 1 1 0 0 0 0 

となります。心の目で六角形グリッドを重ねてください。

一方で、maskを次のように定義します。

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

左上、右下を０にして、それ以外を１としています。Hexマップだと、あるマスに隣接しているのは左右のマスと、上と下の行の２つなので、それを反映しています。
まずは周期境界条件を課して相関を計算してみます。

    N = signal.correlate2d(grid, mask, mode="same", boundary="wrap")

このNは次のようになります。

2 4 2 2 2 1 0 0 0 1 
 3 3 4 1 0 0 1 1 0 2 
  4 3 2 1 1 2 1 2 2 2 
   2 3 2 1 1 2 2 1 3 2 
    2 2 1 1 1 0 2 2 1 1 
     3 4 1 0 1 3 2 1 0 1 
      3 3 2 1 3 5 3 0 0 2 
       4 3 2 3 4 3 1 0 0 1 

自身を含めて周り7マスの中で1の数を数えることができました。
一方でboundaryに"fill"を指定してやると、周りにfill_value（デフォルトで0）を埋めた状態で計算します。つまり周りを壁で囲われている状態、あるいは孤立系です。

    N = signal.correlate2d(grid, mask, mode="same", boundary="fill")

今度のNは次のようになります。

1 3 2 1 0 0 0 0 0 0 
 2 3 4 1 0 0 1 1 0 1 
  2 3 2 1 1 2 1 2 2 2 
   1 3 2 1 1 2 2 1 3 2 
    2 2 1 1 1 0 2 2 1 1 
     3 4 1 0 1 3 2 1 0 0 
      3 3 2 1 3 5 3 0 0 0 
       3 2 1 2 4 3 1 0 0 0 

マインスイーパ

まずはWindows標準で中級の盤面づくり。

    H = 16
    W = 16
    nMine = 40 
    grid = np.zeros(H*W, dtype=int)
    grid[:nMine] = 1
    np.random.shuffle(grid)
    grid = grid.reshape((H, W))

grid行列のなかで１のところに地雷があります。
続いてmaskの定義。今度は自身のマスを含めませんので、中心は0になります。

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

孤立系として１の数を数えます。

    N = signal.correlate2d(grid, mask, mode="same", boundary="fill")

あとは出力ですが、地雷があるマスをmで、数字が０の安全地帯を_で表示します。文字列型のnumpy行列を使います。

    M = np.array(N, dtype=str)
    M[N == 0] = "_"
    M[grid == 1] = "m"
    np.savetxt(sys.stdout, M, "%s")

以下は結果です。まずは地雷マップのgrid。

0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0
1 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0
0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0
1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0
0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0
0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0
0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0
0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0

それに対して盤面の数字。

1 1 _ 1 m 2 2 2 1 _ 1 2 m 1 _ _
m 1 _ 1 1 2 m m 2 _ 1 m 2 1 _ _
2 2 1 _ _ 1 3 m 3 1 1 2 2 1 _ _
1 m 2 1 1 1 2 2 m 2 1 2 m 1 _ _
2 4 m 2 1 m 1 1 2 3 m 2 2 2 1 _
m 3 m 2 1 1 1 _ 1 m 2 1 1 m 1 _
2 3 2 1 1 1 1 _ 1 2 2 1 1 1 1 _
1 m 1 _ 1 m 2 1 _ 1 m 1 1 1 1 _
1 1 1 _ 2 3 m 1 1 2 2 1 1 m 1 _
1 1 1 _ 1 m 3 2 1 m 1 _ 1 2 2 1
1 m 2 1 1 2 m 2 2 2 2 1 1 1 m 1
1 3 m 2 _ 2 2 3 m 2 3 m 2 1 1 1
_ 2 m 2 _ 1 m 2 1 2 m m 3 _ _ _
1 2 1 1 _ 1 1 1 _ 1 4 m 3 _ _ _
m 1 _ _ 1 1 1 _ 1 2 4 m 2 _ _ _
1 1 _ _ 1 m 1 _ 1 m m 2 1 _ _ _

五目並べ

今回はmaskを4種類定義します。

    mask_h = np.zeros((5, 5), dtype=int)
    mask_h[2, :] = 1
    mask_s = np.eye(5, dtype=int)

このmask_hとmask_sの2種類は出力すると次のようになります。

0 0 0 0 0
0 0 0 0 0
1 1 1 1 1
0 0 0 0 0
0 0 0 0 0
----
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1

残り2種類については、mask_hの転置行列と、mask_sの行または列を反転させたものとして定義します。

9x9の盤面について、自分しか石を置いていない状態の盤面を以下のように想定します。

    H = 9
    W = 9
    G = np.zeros((H, W), dtype=int)
    G[0:5, 4:9] = np.eye(5, dtype=int)
    G[4:9, 0:5] = np.eye(5, dtype=int)[:, -1::-1]
    G[1, 4:9] = 2
    G[4:9, 6] = 1
    G[0:3, 0:3] = 1
    np.savetxt(sys.stdout, G, "%d")

1 1 1 0 1 0 0 0 0
1 1 1 0 1 1 1 1 1
1 1 1 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 1 0 1 0 1
0 0 0 1 0 0 1 0 0
0 0 1 0 0 0 1 0 0
0 1 0 0 0 0 1 0 0
1 0 0 0 0 0 1 0 0

1が自分の石のあるところ、0が石のないところ。

これに対して、以下の４種類のmaskに対して、boundary="fill"（つまり孤立系）でcorrelation2dを計算してみます。

    #横一列
    N = signal.correlate2d(G, mask_h, mode="same", boundary="fill")
    np.savetxt(sys.stdout, N, "%d")
    #縦一列
    N = signal.correlate2d(G, mask_h.T, mode="same", boundary="fill")
    np.savetxt(sys.stdout, N, "%d")
    #左上から右下（単位行列）
    N = signal.correlate2d(G, mask_s, mode="same", boundary="fill")
    np.savetxt(sys.stdout, N, "%d")
    #右上から左下
    N = signal.correlate2d(G, mask_s[:, -1::-1], mode="same", boundary="fill")
    np.savetxt(sys.stdout, N, "%d")

これらの結果は以下のようになります。

  mask                 N

----
                3 3 4 3 2 1 1 0 0
                3 3 4 4 4 4 5 4 3
0 0 0 0 0       3 3 3 2 2 1 1 1 1
0 0 0 0 0       0 0 0 0 0 1 1 1 1
1 1 1 1 1       0 0 1 1 2 2 3 2 2
0 0 0 0 0       0 1 1 1 2 2 1 1 1
0 0 0 0 0       1 1 1 1 2 1 1 1 1
                1 1 1 1 1 1 1 1 1
                1 1 1 0 1 1 1 1 1

----

                3 3 3 0 2 1 2 1 1
                3 3 3 0 2 1 2 2 1
0 0 1 0 0       3 3 3 0 3 1 3 2 2
0 0 1 0 0       2 2 2 1 2 1 4 2 2
0 0 1 0 0       1 1 2 1 1 0 4 1 1
0 0 1 0 0       0 1 1 1 1 0 4 1 1
0 0 1 0 0       1 1 1 1 1 0 5 0 1
                1 1 1 1 0 0 4 0 0
                1 1 1 0 0 0 3 0 0

----

                3 2 1 1 3 1 1 1 0
                2 3 2 1 1 4 1 1 1
1 0 0 0 0       1 2 4 2 2 1 5 1 1
0 1 0 0 0       0 2 2 3 2 1 1 4 1
0 0 1 0 0       1 0 2 1 3 1 1 0 3
0 0 0 1 0       0 1 0 1 1 2 1 1 0
0 0 0 0 1       1 0 1 0 2 1 2 1 1
                0 1 0 1 0 2 1 1 1
                1 0 1 0 1 0 1 1 1

----

                1 2 3 2 2 1 1 1 2
                2 3 2 2 1 1 1 2 1
0 0 0 0 1       3 2 2 1 1 1 3 1 2
0 0 0 1 0       2 1 1 1 1 4 1 2 1
0 0 1 0 0       1 0 0 0 4 0 2 1 2
0 1 0 0 0       0 0 0 4 0 2 1 2 1
1 0 0 0 0       0 0 5 0 1 1 2 1 1
                0 4 0 0 1 1 1 1 0
                3 0 0 0 1 1 1 0 0

確かに、maskごとに別の方向について、５つ石が並んだところが5として検出できました。

もし対戦相手も同時にしたいときは、相手の石を-1として、-5となるところを探せば相手が５つ並べたところになるのではないでしょうか。

その他

今回実験に使ったpythonのシェルスクリプトは
https://gist.github.com/sage-git/ad8a6712c96c0eee9f6def2a654ef3a1
にも置いています。

プログラムのパートは以上です。あとはひとりごとです。

数学

僕の記憶では、相関関数は連続・離散関数それぞれで

C(\tau) = \int^\infty_{-\infty}f(t)g(t + \tau)\mathrm{d}t \\
C_t = \frac{1}{T}\sum_{i=0}^T f_i \cdot g_{i + t}

みたいな感じで定義できて、僕は一時期これを使っていろいろな統計力学量を計算していました。一方でConvolution（畳み込み）は

f*g(x) = \int^\infty_{-\infty}f(u)g(x-u)\mathrm{d}u \\
(f*g)_m = \sum_n f_n \cdot g_{m - n} 

と定義されるようです。
相関関数、畳込み共に２次元に対しても容易に拡張が可能です。

特に離散値の定義については、行列の各種操作や画像に対して応用できます。
実際、今回はマス目を数える操作に使ってみました。

画像処理からの発想

別の発想では、画像フィルタの演算も同様の応用ができそうです。
また、流行りのニューラルネットワークについても、畳み込み層も今回のような操作が可能でしょう。
ただの画像のテンプレートマッチだという噂も。

その他

これは決して効率の良い方法ではないでしょう。
数を数えるなら、マス目の数字に特化したアルゴリズムでも考えたほうがいいように思います。

ただpythonのforループを使わなくて良いというのは利点かなと。
ライブラリによっては高度に並列化されてて、うまく使えばGPU化も勝手にできるのではと思います。
でも、なんか鶏を割くに焉んぞ牛刀を用いみがあります。