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`の要素として表現します。

## 配列の出力

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

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

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

乱数を使います。

```fortran
   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$ は初期状態で生存しているセルを生成する確率になります。

これを使って、

```fortran
   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`動かしたものは

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

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

```fortran
   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でライフゲームにおける次の世代の生死条件の定義を確認します。

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

https://ja.wikipedia.org/wiki/ライフゲーム

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

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

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

とできます。

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

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

`where`を使うことで、

```fortran
do i = 1,H
   do j = 1,W
      if( neighbors(i,j) == 3) then
      endif
   enddo
enddo
```
みたいなインデント地獄から逃げました。

# モジュール化プログラム

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

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

```fortran
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
```

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


```fortran
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`が生成されます。

![field.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/16b0efc7-f810-fa26-f858-42ac9cff079d.gif)