ParticleLife, Ising, LifeGame + ChatGPT + DNN

目的

最近TwitterでParticleLifeが話題になっているのを見た。
これって宇宙そのままじゃない？ってすごい刺激を受けた。

Oh, I see. So the branching of the phase transition is too complex and the diversity is expressed by the variation of the initial configuration. https://t.co/e38Q34pK97

— ぼやーじゅ (@voyage_M) March 17, 2023

多体問題に解析解が存在しないっていう事の意味がやっと腑に落ちた #particlelife https://t.co/b5WWe4MTII

— ぼやーじゅ (@voyage_M) March 4, 2023

これのシミュレーション（普通に解析をするのではなく，系をざっくりとらえて、低コストで未来予測する）をやってみたいなと思った。
私たちの世界の複雑さって、この問題にすでに表れているんじゃないかな，と。
あと，こういうのは流行りの深層学習に近い問題なのではと直感した．

でも、取っ掛かりとしてParticleLifeのコードを書いてみたけど，分析も含めると色々難しそうだったので、まずは物理＆多体系と言えば統計力学、Ising模型ってことで、Ising模型のシミュレーションをやってみた。

Ising模型のシミュレーション

ChatGPTにお願いしてみた

ChatGPTにIsing模型のシミュレーションを行うように依頼した．
境界条件やアニメーション出力など，細かくやり取りしたのちの結果で，Google Colabで動作することを確認．

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib import rc
from IPython.display import HTML

def initial_state(N):
    return np.random.choice([-1, 1], size=(N, N))

def energy(spin_config, J, B):
    N = spin_config.shape[0]
    E = -J * np.sum(spin_config * (np.roll(spin_config, 1, axis=0) + np.roll(spin_config, 1, axis=1))) - B * np.sum(spin_config)
    return E

def delta_energy(spin, spin_neighbors, J, B):
    return 2 * J * spin * np.sum(spin_neighbors) + 2 * B * spin

def metropolis(spin_config, J, B, T):
    N = spin_config.shape[0]
    i, j = np.random.randint(N, size=2)
    spin = spin_config[i, j]
    spin_neighbors = [spin_config[(i + 1) % N, j], spin_config[(i - 1) % N, j], spin_config[i, (j + 1) % N], spin_config[i, (j - 1) % N]]
    dE = delta_energy(spin, spin_neighbors, J, B)
    if dE <= 0 or np.random.rand() < np.exp(-dE / T):
        spin *= -1
    spin_config[i, j] = spin
    return spin_config

def animate_simulation(simulation_func, N=50, J=1, B=0, T=1, num_steps=50000, interval=100):
    fig, ax = plt.subplots()
    spin_config = initial_state(N)

    def update(frame):
        nonlocal spin_config
        for _ in range(interval):
            spin_config = simulation_func(spin_config, J, B, T)
        ax.clear()
        ax.imshow(spin_config, cmap="gray")

    ani = animation.FuncAnimation(fig, update, frames=num_steps // interval, interval=interval)
    return ani

if __name__ == "__main__":
    ani = animate_simulation(metropolis)
    rc('animation', html='jshtml')
    ani
    ani.save("gas.gif", writer="pillow")
    #HTML(ani.to_jshtml())

解析結果

一応Ising模型でよく見るような状態が出力た．
ただ，これからこれをDNNでシミュレーションすることを考えているのだけど・・・ちょっと結果が単純すぎて面白味が足りないかな？

以下の話も示唆的で，Ising模型は統計力学の基礎ではあるけど，ParticleLifeみたいなものに現れている不思議さがどこまで表現されるのか疑問にも思った．
ただ，統計力学的な物理量はIsing模型以外を検討するときにも参考になると思った．

非平衡統計力学の研究者としていうと、期待してもらえるのは大変ありがたいけど、正直なところ機械学習に対して大したことは出来ないだろうとかなり厳しい見通しを持っている。一言でいうと、このあたりの非平衡理論は単純化の度合いが高すぎる。https://t.co/mSNbfCUL2i

— シータ (@Perfect_Insider) March 21, 2023

LifeGameのシミュレーション

昔LifeGame（LG）が流行ったことがあったけど，Ising模型よりは難しく，ParticleLifeよりは簡単そうなのでとっかかりにはいいのかなと思った．決定論的なのもとっかかりとしてよい．

ChatGPTに頼んでみた

こっちはなかなか苦戦した．

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

from matplotlib import animation, rc
from IPython.display import HTML

# 100x100のグリッドを作成
grid_size = (100, 100)
grid = np.zeros(grid_size)

# 初期状態を設定
grid[50, 50] = 1
grid[50, 51] = 1
grid[50, 52] = 1

# ライフゲームのルールを定義
def apply_rules(grid):
    # 隣接する8つのセルの状態を取得する
    neighbors = np.zeros(grid_size)
    neighbors[1:, :] += grid[:-1, :]    # 上
    neighbors[:-1, :] += grid[1:, :]    # 下
    neighbors[:, 1:] += grid[:, :-1]    # 左
    neighbors[:, :-1] += grid[:, 1:]    # 右
    neighbors[:-1, :-1] += grid[1:, 1:]    # 左上
    neighbors[:-1, 1:] += grid[1:, :-1]    # 右上
    neighbors[1:, :-1] += grid[:-1, 1:]    # 左下
    neighbors[1:, 1:] += grid[:-1, :-1]    # 右下
    
    # ライフゲームのルールを適用する
    new_grid = np.zeros(grid_size)
    new_grid[(neighbors == 3) | (grid == 1) & (neighbors == 2)] = 1
    
    return new_grid

# アニメーションを作成するための関数
def animate(frame):
    global grid
    plt.clf()
    plt.imshow(grid, cmap='binary')
    plt.axis('off')
    grid = apply_rules(grid)
    return plt

# アニメーションを実行する
fig = plt.figure()
anim = FuncAnimation(plt.gcf(), animate, frames=None, repeat=True, save_count=100, cache_frame_data=False)
rc('animation', html='jshtml')
anim

簡単なDNNで予測をさせてみた

もちろん，ChatGPTにやってもらった．

import os
import glob
import random
from PIL import Image
import numpy as np
from sklearn.model_selection import train_test_split
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Conv2D, Conv2DTranspose, MaxPooling2D

# 予測時間間隔
dt = 3

# 画像を読み込んでNumPy配列に変換する関数
def load_image(path):
    img = Image.open(path).convert('L').resize((100, 100))
    return np.asarray(img) / 255.0

# 入力画像と出力画像のファイルパスのリストを作成する関数
def make_file_list(path, n_samples):
    file_list = os.listdir(path)
    file_list = [os.path.join(path, file) for file in file_list if file.endswith('.png')]
    random.shuffle(file_list)
    return file_list[:n_samples]

# 入力画像と出力画像のデータセットを作成する関数
def make_dataset(file_list, dt):
    x_data = []
    y_data = []
    for i, path in enumerate(file_list):
        img = load_image(path)
        if 5 < i and i + dt < len(file_list):
            x_data.append(img)
            y_data.append(load_image(file_list[i + dt]))
    x_data = np.array(x_data)
    y_data = np.array(y_data)
    return x_data, y_data

# 画像セットのディレクトリパス
path = 'frames2'

# 学習に使用するサンプル数
n_samples = 10000

# ファイルパスのリストを作成する
file_list = []
for seed in range(1, 100):
    dirname = f"{path}/seed_{seed:08d}"
    file_list += make_file_list(dirname, n_samples // 200)

# 入力画像と出力画像のデータセットを作成する
X_train, y_train = make_dataset(file_list, dt)

# データを訓練用と検証用に分割する
X_train, X_valid, y_train, y_valid = train_test_split(X_train, y_train, test_size=0.2, random_state=123)

# モデルの構築
model = Sequential()
model.add(Conv2D(32, (5, 5), activation='relu', padding='same', input_shape=(100, 100, 1)))
model.add(Conv2D(32, (3, 3), activation='relu', padding='same'))
model.add(Conv2DTranspose(1, (3, 3), activation='sigmoid', padding='same'))
model.compile(loss='binary_crossentropy', optimizer='adam')

# モデルの学習
model.fit(X_train, y_train, epochs=50, batch_size=32, validation_data=(X_valid, y_valid))

# モデルの保存
model.save('my_model.h5')

解析結果

学習データセットとして，初期条件をランダムに変えたLGのシミュレーションを行い，ある時刻tの画像とdt（＝3）後の画像の組を作成する．
前者を入力，後者を出力として適当なDNNを構成して学習させる．
その後，DNNで予測した白黒の分布を正解と比較する．
予測結果はfloat値で得られるので，0.5を境に2値化したものも併せて示す．

ちょっとやってみた割には，思ったよりそれらしい結果が得られたと思う．

• 周辺に何もない固定物体については比較的正しく予想できている．

• 十字と円形を周期的に移り変わる物体についても，きちんと予測できているものもありそう．

• 最初に大き目の畳み込み層を設定しないと，そのサイズを超えるような構造は抽出できないだろう．

• 最初に設定する畳み込み層を大きくするほど予測の精度は上がるが学習も大変そう．

今後の展望

DNNの層の構成を変えて検討してみたい．
例えば，畳み込み層としてN×Mの層を設定すれば，そのサイズの構造を抽出できる．
初期条件としてN×Mしか値を持たなかったらそれだけで時間発展は追えるだろう．
実際はN×Mの横にも隣接するN’×M’のセルが存在し，それとの相互作用で時間発展は決まる．
という事は，サイズの異なる複数の層を設定し，それを全結合したほうが予測精度は上がるだろうか・・・など，DNNの構成を変えて予測がどう変わるか見てみたい．
また，dtを変更して結果がどう変わるかなども気になるところである．

最後に

というわけで，今回の検討は以上です．
時間を見てもう少し検討してみたいと思います．
バグとか面白い案があったら教えてください．