9-9. vcoptで逆ライフゲーム | Vignette & Clarity（ビネット＆クラリティ）

9-9. vcoptで逆ライフゲーム

vcopt応用

2019/06/23

やること
参考文献
実行環境
ライフゲームの基本
ソースコード
結果
順方向で確認

やること

ライフゲームは１ステップ先を計算するのは簡単ですが、１ステップ前を計算するのはなかなか難しいです。vcoptで逆ライフゲームを試してみましょう。今回はvcopt().dcGA()を繰り返し実行するのでトリッキーです。

参考文献

ライフゲームは、順方向の結果は一意に定まりますが、逆方向は一意ではありません。したがって、ハッシュや暗号に応用できると言われています。

https://www.jstage.jst.go.jp/article/pjsai/JSAI08/0/JSAI08_0_196/_pdf

vcoptの使い方についてはチュートリアルをご参照ください。

vcoptの仕様については最新の仕様書をご参照ください。本記事執筆時とは仕様が異なる場合があります。

実行環境

WinPython3.6をおすすめしています。

WinPython - Browse /WinPython_3.6/3.6.7.0 at SourceForge.net

Portable Scientific Python 2/3 32/64bit Distribution for Windows

ライフゲームの基本

こちらの記事をご参照ください。

ソースコード

長い。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from vcopt import vcopt

#============================
#設定
#============================

#高さ、幅
h, w = 10, 10

#任意の最終状態を用意
state = np.array([[0,0,1],
                  [1,0,1],
                  [0,1,1]])

#フィールドのどこに置くか（左上点を指定）
p = (4, 4)











#============================
#関数など
#============================
#目標フィールドの生成
f_aim = np.zeros((h, w), dtype=bool)

#任意の状態を置く
f_aim[p[0]:p[0]+len(state), p[1]:p[1]+len(state[0])] = state

#畳み込み用のフィルタ
g = np.array([[1.0, 1.0, 1.0],
              [1.0, 0.0, 1.0],
              [1.0, 1.0, 1.0]])

h_crop, w_crop = None, None
p_crop = (None, None)

#paraからfを復元する関数
def para_to_f(para, h_crop, w_crop, p_crop):
    #クロップ部分を復元
    f_crop = para.reshape(h_crop, w_crop)
    #フィールド全体の復元
    f = np.zeros((h, w), dtype=bool)
    f[p_crop[0]:p_crop[0]+h_crop, p_crop[1]:p_crop[1]+w_crop] = f_crop
    
    return f

#fからfutureを計算する関数
def f_to_future(f):
    #周囲の生存マスを足し込む
    mask = signal.convolve2d(f, g, mode='same', boundary='wrap')
    
    #１ステップ未来のフィールド（すべて死状態）
    future = np.zeros(f.shape, dtype=int)
    
    #生きているマスが生きる条件（＝生存）
    future[mask*f==2] = 1
    future[mask*f==3] = 1
    #死んでいるマスが生きる条件（＝誕生）
    future[mask*~f==3] = 1
    
    return future

#３つのグラフを並べて表示する関数
def show_3(f1, f2, f3, title1, title2, title3, i):
    plt.figure(figsize=(10, 10))
    
    plt.subplot(1, 3, 1)
    plt.imshow(f1, cmap='inferno')
    plt.title(title1)
    plt.subplot(1, 3, 2)
    plt.imshow(f2, cmap='inferno')
    plt.title(title2)
    plt.subplot(1, 3, 3)
    plt.imshow(f3, cmap='inferno')
    plt.title(title3)
    plt.savefig('save/{}.png'.format(i), bbox_inches='tight', pad_inches=0)
    plt.show(), print()

#生きているセルのある部分を余白つけてクロップする関数
def crop(f):
    col_sum = np.sum(f, axis=0)
    row_sum = np.sum(f, axis=1)
    
    col_exist = np.where(col_sum>0)[0]
    col_start = max(col_exist[0] - 1, 0)
    col_end = min(col_exist[-1] + 2, f.shape[1])
    
    row_exist = np.where(row_sum>0)[0]
    row_start = max(row_exist[0] - 1, 0)
    row_end = min(row_exist[-1] + 2, f.shape[0])
    
    f_crop = f[row_start:row_end, col_start:col_end]
    
    h_crop, w_crop = f_crop.shape[0], f_crop.shape[1]
    p_crop = (row_start, col_start)
    
    return f_crop, h_crop, w_crop, p_crop

#評価関数
def rGL_score(para):
    
    #paraからfの復元
    f = para_to_f(para, h_crop, w_crop, p_crop)
    
    #１ステップ未来のフィールド
    future = f_to_future(f)
    
    #ステップ前後の差分マス数
    diff = np.sum((future - f_aim)**2)
    
    #前ステップの生セル比率
    l_score = np.sum(f) / (h * w)
    
    #スコアを返す
    return diff + l_score


#============================
#リバース開始
#============================

score = 0.0
step = -1
while score < 1.0:
    
    print('{}stepを計算中'.format(step))

    #生きているセルのある部分を余白つけてクロップ
    f_crop, h_crop, w_crop, p_crop = crop(f_aim)

    #パラメータ範囲
    para_range = [[i for i in range(2)] for j in range(h_crop * w_crop)]
    
    #GAで最適化
    para, score = vcopt().dcGA(para_range,
                                rGL_score,
                                0.0,
                                show_pool_func='bar',
                                seed=1,
                                pool_num=1000)
    #結果の表示
    print('score={}'.format(score))
    
    #paraからfの復元
    f = para_to_f(para, h_crop, w_crop, p_crop)
    
    #１ステップ未来のフィールド
    future = f_to_future(f)
    
    #表示
    show_3(f, future, f_aim, 'before', 'after', 'f_aim', -step)

    #更新
    f_aim = f
    step -= 1

評価関数を少し工夫しています。あるparaからフィールドを復元し、１ステップ後を計算し、目標フィールドとの差分をdiffとします。これだけだとセルが拡散しがちなので、セルがまとまりやすいようにゴニョゴニョしています（説明が面倒くさくなり絶命）。

結果

ソースコードをそのまま実行すると、グライダーの逆ライフゲーム探索が始まります。

-1stepを計算中
 [|<                                                ] gen=73000, best_score=0.05
score=0.05

-2stepを計算中
 [  |<                                              ] gen=55000, best_score=0.05
score=0.05

-3stepを計算中
 [|<                                                ] gen=87000, best_score=0.05
score=0.05

-4stepを計算中
 [|<                                                ] gen=77000, best_score=0.06
score=0.06

-5stepを計算中
 [|<                                                ] gen=71000, best_score=0.08
score=0.08

-6stepを計算中
 [|<                                                ] gen=143000, best_score=0.1
score=0.1

-7stepを計算中
 [|<                                                ] gen=138000, best_score=0.33
score=0.33

-8stepを計算中
 [    |<                                            ] gen=194000, best_score=3.23
score=3.23

ということで、-7ステップ目で急激に発散し、-8ステップ目は見つかりませんでした。score=3.23というのは、３マスの食い違いが生じ、フィールドの23%が生セルという意味です。

順方向で確認

-7ステップ目から再生してみます。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

#============================
#初期状態の設定
#============================

#高さ、幅
h, w = 10, 10

#任意の状態を用意
state = np.array([[0,0,1,0,1,1,1],
                  [0,1,1,0,1,0,0],
                  [1,0,1,1,0,0,0],
                  [0,0,1,1,1,1,1],
                  [0,1,0,0,0,1,0],
                  [0,0,1,1,0,1,0],
                  [0,1,1,0,1,0,0],
                  [0,0,1,0,1,1,0],
                  [1,0,1,1,0,1,0],
                  [0,0,1,1,1,0,0]])

#フィールドのどこに置くか（左上点を指定）
p = (0, 1)

#終了ステップ数
max_step = 15

#============================
#メイン処理
#============================

#フィールドの生成
f = np.zeros((h, w), dtype=bool)

#任意の状態を置く
f[p[0]:p[0]+len(state), p[1]:p[1]+len(state[0])] = state

#畳み込み用のフィルタ
g = np.array([[1.0, 1.0, 1.0],
              [1.0, 0.0, 1.0],
              [1.0, 1.0, 1.0]])

#初期状態の表示
#plt.figure(figsize=(10, 10))
plt.imshow(f, cmap='inferno')
#plt.savefig('save/{}.png'.format(0), bbox_inches='tight', pad_inches=0)
plt.show(), print()

#状態の更新
for i in range(1, max_step + 1):
    
    #周囲の生存マスを足し込む
    mask = signal.convolve2d(f, g, mode='same', boundary='wrap')
    
    #１ステップ未来のフィールド（すべて死状態）
    future = np.zeros(f.shape, dtype=bool)
    
    #生きているマスが生きる条件（＝生存）
    future[mask*f==2] = 1
    future[mask*f==3] = 1
    #死んでいるマスが生きる条件（＝誕生）
    future[mask*~f==3] = 1
    
    #フィールドの更新（浅いコピーに注意）
    f = future
    
    #表示
    #plt.figure(figsize=(10, 10))
    plt.imshow(f, cmap='inferno')
    #plt.savefig('save/{}.png'.format(i), bbox_inches='tight', pad_inches=0)
    plt.show(), print()

結果