Fortranでじゃんけん型セル・オートマトン

これまでセル・オートマトンで遊ぶ記事をいくつか書いていたのですが、カレンダーの穴埋めに、この流れでもうひとつ。

「 Chainerとcupyをcuda版numpyとして使ってオートマトン」のモデルを「Fortranでライフゲーム」を踏まえてFortranで実装します。

セル・オートマトン

セル・オートマトンは、以下を定義したセルを並べて時間発展させることによって遊ぶことのできるモデルです。

• 個々の状態
• 他のセルとの相互作用で自分の状態がどう変わるか

プログラムではこのルールしか決めていないのに、集団になると多彩な現象をみせる不思議さが人気の秘密です。

詳しくはWikipediaなどをご参照ください。
Wikipediaでは、状態が生死の２種で23/3ルールの俗に言うライフゲーム（Conway's Game of Life）が妙に詳しいです。ライフゲーム中に出てくるパターンそれぞれに専用ページがあるレベル。

ここでは個人的にじゃんけんオートマトンと認識している、３すくみのオートマトンを実装します。
例えば https://qiita.com/SouTakenaka/items/a7dcfbde0a08063f4d41 にあるように、BZ反応のシミュレーションに用いられたりします。

ですが今回は
https://www.youtube.com/watch?v=zhxSyUL3OOY
にあるように、渦を作って遊ぶだけにしておきます。

セルの種類と相互作用

正方格子にセルを並べ、隣接セルをマインスイーパ式（つまり上下左右と斜め４つの系）とします。

それぞれのセルについて以下の４種類を定義します。

また、セルはHPというパラメータを持っています。このパラメータは隣接セルとの相互作用により減らされます。

HPが０になると、セルの状態は上図の矢印の方向に状態が変異します。

なお、石以外の紙・ハサミと隣接したEmptyセルも石になります。複数種類のセルと隣接したEmptyのルールとか、あと「emptyセルがどのセルと隣接していたか」を管理する分プログラムが面倒になりますので、一律に石に変化するとしました。

ちなみにEMPTYは初期配置にしか現れないということにしました。

HPについて

Emptyは他の石・紙・ハサミセルから、また石・紙・ハサミセルはじゃんけんで負けてしまう相手から１つのセルにつき１減らされます。
例えば紙セル４つと隣接した石セルは１ターンでHPを４減らされます。また、ハサミセルを包囲しているすべての紙セルは１ずつHPを減らされます。

変異した直後はHPを１０とします。

初期配置

０〜３の整数をランダムで作って並べます。

実装

アルゴリズムの各段階のスニペットを解説します。全部まとめてのプログラムもあります。

盤面準備

        W = width
        H = height
        allocate(RPS(W, H))
        allocate(HP(W, H))
        RPS(:, :) = 0
        HP(:, :) = 1

RPSがセルの種類を保持します。
HPがセルのHPを保持します。

とりあえず全部のセルをHP1のEMPTYにしています。

初期配置

何の変哲もないループです。

        integer :: x, y
        double precision :: p
        do x = 1, W
            do y = 1, H
                call random_number(p)
                RPS(x, y) = int(p*4.d0)
            enddo
        enddo

数え上げ

「Fortranでライフゲーム」では、周囲のセルをカウントするとき、

        neighbors = 0
        do concurrent (dx = -1:1, dy = -1:1)
            neighbors = neighbors + cshift(cshift(field, dx, 1), dy, 2)
        enddo

としましたが、これは一時的にとはいえcshiftした結果の配列が毎回生成される気がしたので、愚直にループを回すことにしました。
あと、cshiftを用いると周期境界条件が簡単に実現できますが、自作doループ数え上げについては配列端での処理が面倒になります。そのため、今回は端で切れている孤立系です。

        do concurrent (x=1:W, y=1:H)
            do concurrent (dx = max(1,x-1):min(W,x+1), dy = max(1,y-1):min(H,y+1))
                if( dx == x .and. dy == y) cycle
                i = RPS(x, y) + 1
                if(i == 4) i = 1
                if(RPS(dx, dy) == i) neighbor(x, y) = neighbor(x, y) + 1
                if(RPS(x, y) == 0 .and. RPS(dx, dy) > 0) neighbor(x, y) = neighbor(x, y) + 1
            enddo
        enddo

do concurrentは２次元のdoループを１つにまとめるためだけで、特に並列化とかは狙っていません。
一番外側のループはx, yについて回すもので、システム全体をスキャンします。
その内側のdx、dyは注目しているx、yについて隣接するセルをスキャンします。上下左右に±1してループするので９回ループですが、dx = dy = 0となるのは自分自身を参照しているのでcycleでとっとと次に送ります。配列の端に来たときに１ずらすと配列の外を参照してしまう可能性がありますので、minやmaxを使ってはみ出ないようにします。　

iは(x, y)のセルに勝てる種類を整数ラベルで持っています。(dx, dy)がiであれば(x, y)に1ダメージを与えます。
　つまり各セルごとに問題となる種類のセルを数え上げることになるので、これをneighborにたくわえていきます。
　あと、EMPTYセルはEMPTYでないセルから無条件にダメージを喰らいますので、これも別に数えます。

変数のスコープについて。

• x y dx dy i jはループ毎に独立した変数
• RPSはモジュール変数で共通で変更なし
• neighborはループ間で共用ですが集計のみの変更

do concurrentの使い方はこれでいいのだろうか。

セルの変異

周りのセルの種類を数え上げたのちは、それに基づいてHP計算をします。

       HP(:, :) = HP(:, :) - neighbor(:, :)

HPが０になったセルを変異させます。
単純に種類変数を１ふやし、４になったものを１に戻してやれば、３すくみが再現できます。

        logical :: mute(W, H)
        integer, parameter :: MAX_HP = 10
        mute(:, :) = HP(:, :) < 1
        where(mute) RPS(:, :) = RPS(:, :) + 1
        where(RPS == 4) RPS = 1
        where(mute) HP = MAX_HP

画像生成

最後に可視化のための関数。

    subroutine set_image(img)
        character(len=1), intent(inout) :: img(3, W, H)
        integer :: x, y, s
        do concurrent(x=1:W, y=1:H)
            s = RPS(x, y)
            
            if(s == 0) img(:, x, y) = EMPTY_COLOR
            if(s == 1) img(:, x, y) = ROCK_COLOR
            if(s == 2) img(:, x, y) = PAPER_COLOR
            if(s == 3) img(:, x, y) = SCISSOR_COLOR
        enddo
    end subroutine set_image

ここで、色は以下のようにセットしました。

    character(len=1), parameter :: EMPTY_COLOR(3) = char([0, 0, 0])  !B,G,R
    character(len=1), parameter :: ROCK_COLOR(3) = char([12, 0, 179])  !B,G,R
    character(len=1), parameter :: PAPER_COLOR(3) = char([216, 216, 216])
    character(len=1), parameter :: SCISSOR_COLOR(3) = char([0, 179, 44])

クリスマスに合わせて、赤・緑・白のクリスマスカラーにしました。

• EMPTYは黒。
• ROCKは赤。サンタの服とかトナカイの鼻とか。
• PAPERは白。若干灰というか銀色。雪ですね。
• SCISSORは緑。もみの木の葉のイメージ。

実際の色は「セルの種類と相互作用」節のグラフをご参照ください。ノードの背景色を対応する色にしています。

以上がautomatonモジュールにまとめておく処理群になります。

初期化と時間発展

要はメイン関数に書くもの。

    H = 256
    W = 256
    N = 800
    call init_system(H, W) ! 盤面準備
    call random_init()     ! 初期配置
    do i = 1, N
        call step()        ! 数え上げと変異
        ! ここで可視化
    enddo

このように、とりあえず256x256の盤面で800ステップ時間発展させてみます。

ソースコード全体

以上のアルゴリズムを、Fortranのプログラムとして実行できるように整え、またOpenCVを使ってリアルタイムで表示できるようにします。

automaton.f90

module rps_automaton
    implicit none
    integer, allocatable :: RPS(:, :), HP(:, :)
    integer :: W, H
    integer, parameter :: MAX_HP = 10
    character(len=1), parameter :: EMPTY_COLOR(3) = char([0, 0, 0])  !B,G,R
    character(len=1), parameter :: ROCK_COLOR(3) = char([12, 0, 179])  !B,G,R
    character(len=1), parameter :: PAPER_COLOR(3) = char([216, 216, 216])
    character(len=1), parameter :: SCISSOR_COLOR(3) = char([0, 179, 44])

contains
    subroutine init_system(height, width)
        integer, intent(in) :: height, width
        W = width
        H = height
        allocate(RPS(W, H))
        allocate(HP(W, H))
        RPS(:, :) = 0
        HP(:, :) = 1
    end subroutine init_system

    subroutine step()
        integer :: neighbor(W, H)
        integer :: x, y, dx, dy, i, j
        neighbor(:, :) = 0
        do concurrent (x=1:W, y=1:H)
            do concurrent (dx = max(1,x-1):min(W,x+1), dy = max(1,y-1):min(H,y+1))
                if( dx == x .and. dy == y) cycle
                i = RPS(x, y) + 1
                if(i == 4) i = 1
                if(RPS(dx, dy) == i) neighbor(x, y) = neighbor(x, y) + 1
                if(RPS(x, y) == 0 .and. RPS(dx, dy) > 0) neighbor(x, y) = neighbor(x, y) + 1
            enddo
        enddo
        HP(:, :) = HP(:, :) - neighbor(:, :)
        call mutation(HP < 1)
    end subroutine step

    subroutine mutation(mute)
        logical, intent(in) :: mute(W, H)
        where(mute) RPS(:, :) = RPS(:, :) + 1
        where(RPS == 4) RPS = 1
        where(mute) HP = MAX_HP
    end subroutine mutation

    subroutine random_init()
        integer :: x, y
        double precision :: p
        do x = 1, W
            do y = 1, H
                call random_number(p)
                RPS(x, y) = int(p*4.d0)
            enddo
        enddo
    end subroutine random_init

    subroutine set_image(img)
        character(len=1), intent(inout) :: img(3, W, H)
        integer :: x, y, s
        do concurrent(x=1:W, y=1:H)
            s = RPS(x, y)
            
            if(s == 0) img(:, x, y) = EMPTY_COLOR
            if(s == 1) img(:, x, y) = ROCK_COLOR
            if(s == 2) img(:, x, y) = PAPER_COLOR
            if(s == 3) img(:, x, y) = SCISSOR_COLOR
        enddo
    end subroutine set_image

end module rps_automaton

main.f90

program main
    use rps_automaton, only: init_system, step, random_init, set_image
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    use iso_c_binding
    implicit none
    integer :: N, i
    integer(c_int) :: H, W, ret
    character(len=1), pointer :: img(:, :, :)
    type(c_ptr) :: p_img
    type(CvMat), pointer :: cvimg

    H = 256
    W = 256
    N = 800

    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_8U, 3))
    call c_f_pointer(p_img, cvimg)
    call c_f_pointer(cvimg%data, img, (/3, W, H/))

    call init_system(H, W)
    call random_init()
    do i = 1, N
        call step()
        call set_image(img)
        call cvShowImage("hoge", p_img)       
        ret = cvWaitKey(10)
    enddo
    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)

end program main

あと、FortranからOpenCVを叩く　画像処理関数編を参考¹に、core_c.f90とhighgui_c.f90、types_c.f90を同じフォルダに準備します。

コンパイルと実行

gfortran types_c.f90 core_c.f90 imgcodecs_c.f90 highgui_c.f90 automaton.f90 main.f90 -L/usr/local/lib -lopencv_core -lopencv_highgui
./a.out

Fortranのソースコードを一度にコンパイルするときは、引数でのファイルの順番に気をつけてください。例えばtypes_c.f90をcore_c.f90の後に持ってくると、「Fatal Error: Can't open module file ‘cv_types_c.mod’ for reading at (1): そのようなファイルやディレクトリはありません」と言われます。
OpenCVのGUIウィンドウが立ち上がり、セル・オートマトンの様子がアニメーションで見れるはずです。

スナップショット

以下の画像を用意するにあたり、デスクトップのスクリーンキャプチャではうまく行かなかったので、
FortranからOpenCVを叩く　画像処理関数編のcvImWriteを使いました。

• 1フレーム目
  

• 6フレーム目
  

• 12フレーム目
  

• 20フレーム目
  

• 30フレーム目
  

• 45フレーム目
  

• 70フレーム目
  

• 120フレーム目
  

• 500フレーム目
  

• 800フレームまでのアニメーション
  

ファイルサイズの関係で、フレームを間引きしています。

補足

このセルと相互作用の定義に行き着くまでに、六角格子ですが前にQiitaに投稿した「人工生命を実装して世界大戦、人類滅亡への一歩を踏み出す」で試行錯誤しています。特にHPシステムとか。

ちなみにこの六角格子とかcupy版とかですが、セルの突然変異とか入れているので、決定論的ではないです。初期配置でセル３種類の配置で勝負し、突然変異なしというのが本来のセル・オートマトンかなと思います。まあ突然変異確率を減衰するようにしておりますので、これが十分小さくなるまではセル・オートマトンの初期配置生成段階、あるいは撹拌を再現していると思えばいいかなと。

その他

cupy版のやつで大きいサイズで長時間動かしてたら、渦の中心がゆっくり移動しているように見えたのですが、これはどんなモデルで表現できるんだろうか。

終わりに

クリスマスらしい色合いにしましたので、クリスマスらしさを感じられていただけたら幸いです。

1. highgui_c.f90の実装において、cvWaitKey(msec)のmsecにvalue属性を付け忘れていました。すみません。 ↩