Fortranでライフゲーム

この前のAPLの記事の参考で書いた、APLのライフゲームの動画。
https://www.youtube.com/watch?v=a9xAKttWgP4
これを参考にしながら、Fortranでも実装してみます。

このとき、できるだけロジックの行数を短くしようかと思います。

つまり目的はライフゲームそのものよりも、Fortran文法の確認です。

コンパイラ

いろいろやってみたらFortran2003以降の文法も使いましたので、バージョンを書いておきます。

% gfortran --version
GNU Fortran (Ubuntu 5.4.0-6ubuntu1~16.04.9) 5.4.0 20160609

下ごしらえ

まずは必要なサブルーチンを用意します。

今回はセルが並んだ環境をW x Hの整数型の２次元配列fieldで表し、生存しているセルを1の要素、死亡しているセルを0の要素として表現します。

配列の出力

   do i = 1, H
      write(*, '(*(i0,:," "))') field(:, i)
   enddo

FortranのWrite文のformatを少し工夫します。
(*(i0,:," "))は(i0,:," ")が*個あることを示します。
:はデータがなくなれば処理を打ち切り、" "は区切り文字を空白にします。
https://qiita.com/chirimen/items/f867855f8707ed4e104c

どうせ0か1しか出力しないので桁数は勝手に揃います。

盤面作成、初期状態の作成

乱数を使います。

   double precision, allocatable :: s(:, :)
   allocate( s(W, H) )
   call random_number(s)

random_numberサブルーチンにより、配列sにインプレースで$(0, 1)$の乱数が代入されます。

ここで浮動小数点型を整数型に代入すると切り捨てられることを利用します。
$0 < s < 1$ なる乱数 $s$に$0<p<1$の値を足して切り捨てると、$0 < s < 1-p$のときは０に、$ 1-p \leq s < 1 $のときは１になります。
つまり $p$ は初期状態で生存しているセルを生成する確率になります。

これを使って、

   integer, allocatable :: field(:, :)
   double precision :: start_ratio
   allocate( field(W, H) )
   field = s + start_ratio

とすればfieldの中身は

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

のようになります。（先に書きました方法で出力しました。）

周囲の生存数カウント

これはある要素について、周囲８マスの１の個数を数えることに相当します。
ただ冒頭の動画でもそうですが、ここでは自身も含めて数えることとします。

方向性としては上下左右斜め８方向にずらしたものをすべて足し合わせます。
Fortranではcshift(array, shift, dimension)を用います。

cshift関数で動かした様子です。

field
1 1 1 1 0
0 0 0 0 0
0 0 1 0 0
1 0 1 1 0
 ------------------
 cshift(field, 1, 1) (左に１つ動かす)
1 1 1 0 1
0 0 0 0 0
0 1 0 0 0
0 1 1 0 1
 ------------------
 cshift(field, 1, 2) (上に１つ動かす)
0 0 0 0 0
0 0 1 0 0
1 0 1 1 0
1 1 1 1 0

横方向にdx、縦方向にdy動かしたものは

cshift( cshift(field, dx, 1), dy, 2 )

となります。
そしてdx、dyをそれぞれ-1, 0, 1の３通り、合わせて９通り動かしたものを合計します。

   integer, allocatable :: neighbors(:, :) 
   allocate( neighbors(W, H) )
   neighbors = 0
 
   do concurrent (dx = -1:1, dy = -1:1)
      neighbors = neighbors + cshift(cshift(field, dx, 1), dy, 2)
   enddo

concurrentにより、dxとdyの２重ループを１つにまとめました。
concurrentはループの各段階で独立であることが求められていますが、整数の合計なら大丈夫でしょう。
このneighborsはこのようになります。

 field
1 1 1 1 0
0 0 0 0 0
0 0 1 0 0
1 0 1 1 0
------------------
 neighbors
3 5 5 4 4
2 4 4 3 2
1 3 3 3 2
3 6 6 5 4

cshiftが循環的に動かしますので、いきおい今回の実装では周期境界条件となります。

Fortranでなければ、畳み込みもしくはBoxFilterで実装できそうです。
Numpyであればconvolveもしくはcorrelate2dが使えそうです。
https://ja.stackoverflow.com/questions/22363/pythonのconvolve2dについて

次の世代

ここでWikipediaでライフゲームにおける次の世代の生死条件の定義を確認します。

誕生
　　死んでいるセルに隣接する生きたセルがちょうど3つあれば、次の世代が誕生する。
生存
　　生きているセルに隣接する生きたセルが2つか3つならば、次の世代でも生存する。
過疎
　　生きているセルに隣接する生きたセルが1つ以下ならば、過疎により死滅する。
過密
　　生きているセルに隣接する生きたセルが4つ以上ならば、過密により死滅する。

ここで、自身の生存も足したneighbors配列においては、

• 3のとき
  • 自身が死亡なら周囲に生存セルが３つあるので、次の世代では誕生
  • 自身が生存なら周囲に生存セルが２つあるので、次の世代では生存
• ４のとき
  • 自身が死亡なら周囲に生存セルが４つあるので、次の世代では死亡のまま
  • 自身が生存なら周囲に生存セルが３つあるので、次の世代では生存のまま
• それ以外のとき(1、2、5〜9)
  • 自身の生存・死亡にかかわらず、次の世代では死亡

なので、

• neighborsが３の要素は必ず1
• neighborsが４の要素はfieldの値をそのまま継続
• それ以外の要素では必ず0

とできます。

これをFortranコードにすると以下のようになります。

   where(neighbors == 3)
      field = 1
   elsewhere(neighbors == 4)
      field = field
   elsewhere
      field = 0
   endwhere

whereを使うことで、

do i = 1,H
   do j = 1,W
      if( neighbors(i,j) == 3) then
      endif
   enddo
enddo

みたいなインデント地獄から逃げました。

モジュール化プログラム

Fortran2003ではクラスも定義できるのですが、これくらいのスケールならモジュールが一番書きやすいように思います。
ライブラリとして再利用するのならクラス化したほうがいいでしょう。

まずfieldの状態と各種サブルーチンを持つモジュールを作ります。

module lifegame
    implicit none
    integer :: W, H
    integer, allocatable :: field(:, :)
    integer, allocatable :: neighbors(:, :)
contains
    subroutine init_field(width, height, start_ratio)
        integer, intent(in) :: width, height
        double precision, intent(in) :: start_ratio
        double precision, allocatable :: s(:, :)
        W = width
        H = height
        allocate( field(W, H) )
        allocate( s(W, H) )
        allocate( neighbors(W, H) )
        call random_number(s)
        field = s + start_ratio
    end subroutine init_field

    subroutine print_field()
        implicit none
        integer :: n, i
        do i = 1, H
            write(*, '(*(i0,:," "))') field(:, i)
        enddo
    end subroutine print_field

    subroutine one_epoch
        implicit none
        integer :: dx, dy, i

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

        where(neighbors == 3)
            field = 1
        elsewhere(neighbors == 4)
            field = field
        elsewhere
            field = 0
        endwhere
    end subroutine one_epoch
end module lifegame

そして、この下にこのモジュールを使うプログラムを書いていきます。

program flg
    use lifegame
    implicit none
    integer :: epoch, width, height
    integer :: i

    write(*, *) "# Width?"
    read(*, *) width
    write(*, *) "# Height?"
    read(*, *) height
    write(*, *) "# Epoch?"
    read(*, *) epoch

    call init_field(width, height, 0.5d0)
    call print_field

    do i = 1, epoch
        call one_epoch
        write(*, *)
        write(*, *)
        call print_field
    enddo 

    stop
end program flg

これらのコードをまとめてlifegame.f90として保存します。
同じソースコードファイルに書く関係上、上からモジュール、プログラムの順で書く必要があります。

これを実行すると、標準出力に

 # Width?
4
 # Height?
4
 # Epoch?
12  
1 1 1 1
0 0 0 0
0 0 0 0
1 0 0 1


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


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

・・・

という感じで、盤面、空行２行、盤面、空行２行、・・・と表示されます。

可視化

gnuplotを使います。

reset
set key below

set nozeroaxis
set notics
unset colorbox

set size square

set term gif animate size 240,200
set output "field.gif"

set palette grey
do for [n=0:200:1] {
   plot "field.log" index n matrix with image title sprintf("generation %d", n)
}

最近のgnuplotらしくfor構文を導入しました。

これをplotter.pltファイルとして保存します。

世代ごとに２行の空行を入れることで、plotのときにindexを指定して世代を選んでプロットできるようになります。

なおプロット行は

   p "field.log" i n mat w ima t sprintf("generation %d", n)

とキーワードを省略できます。

実行

端末上で次のように実行します。

$ gfortran lifegame.f90
$ echo "64\n64\n512" | ./a.out > field.log
$ gnuplot plotter.plt

最後のgnuplotにより、field.gifが生成されます。