JuliaGPU / ArrayFire でイジングモデルを計算する

# 概要

## 動機

- 実装を通してイジングモデルを理解する
- GPU における計算手法を確立する

GPU 計算については、最高のパフォーマンスを求めるというよりも、できる限り楽に、そこそこのパフォーマンスを出せる方法を模索したい。

## イジングモデルについて

特に参考にさせてもらったサイト

[Pythonで粒子のスピンを操ろう！: 磁石とはなんぞや？](https://qiita.com/jabberwocky0139/items/a34c9e9a2e711d4d6a54)
[Ising Modelを平易に解説してみる](http://enakai00.hatenablog.com/entry/20150106/1420538321)
[Gen Kuroki さんの実装](http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7)

## 結果

[![イジングモデル動画](https://qiita-image-store.s3.amazonaws.com/0/202534/6654033c-d753-e457-b93a-1211c0b3922c.png)](https://imgur.com/8yqpnzt)

画像クリックか下のアドレスでGIF動画（12MB）へ
https://imgur.com/8yqpnzt

コード （Julia Notebook）
https://gist.github.com/Lirimy/127d6a557ba10d549d17c6687ee9c5c2

１フレームですべてのピクセルを平均１回評価する。
サイズ 640×360, 450 フレーム (30 fps * 15 秒)
計算だけなら 3 秒、動画を出力しても 11 秒で完了する。

| サイズ        |       640x360 |      1920x1080 |
|:-------------|--------------:|---------------:|
| 計算のみ      |       2.8 sec |        20 sec |
| 計算＋画像出力 |       7.0 sec |        51 sec |
| 動画化 (gif)  |       3.8 sec |        34 sec |
| 動画化 (mp4)  |       5.8 sec |        52 sec |
| 動画容量 (gif) |        12 MB |        103 MB |
| 動画容量 (mp4) |        33 MB |        298 MB |
| GPU 使用率     |         30 % |          50 % |

<details><summary>イジングモデル動画化についての知見</summary><div>

当初は　mp4 形式で YouTube にアップロードしようと目論んで、 [実際やった](https://www.youtube.com/watch?v=3Pio46Zhgsg) のだが、再エンコードされてしまって残念な画質になってしまっている。アップロード前の元動画は 15秒 33 MB でおよそ 18 Mbps を必要とするが、 [YouTube の推奨ビットレートが 1 Mbps](https://support.google.com/youtube/answer/1722171?hl=ja) なので、当然の成り行きである。

そもそもイジングモデルではピクセルがランダムにフリップするので、空間的に相関が少なく、いわゆるエンコ殺しの画像になっており、mp4とは相性が悪いようだ。一方で、色的に見ると２状態に対応する２色しか使われていないので、gifのパレットが有効に機能するようだ。gif動画はmp4と同条件で 12 MB であり、 およそ　1/3 の容量に収まっている。動画変換も速い。

</div></details>

# 実装の準備

## GPU ライブラリの選択

可視化に優れた Julia を使いたいので、今回は ArrayFire のラッパー ArrayFire.jl を選択した。

- 高レベルなので書きやすい
- バックエンドを CUDA と OpenCL から選べる

今回は見送ったが、他の候補たち

- ArrayFire (C++) を Julia から呼び出し
- CUDA ラッパー
- CUDA Fortran を Julia から呼び出し

## ArrayFire.jl の使い方

ライブラリを使う場合、普通は公式のドキュメントを読めばなんとかなるのだが、残念ながらまだドキュメントが整っていないようなのでどうすべきか、という話。

1. [ArrayFire.jl](https://github.com/JuliaComputing/ArrayFire.jl) の Supported Functions で欲しい関数が実装されているか探す
2. [wrap.jl](https://github.com/JuliaComputing/ArrayFire.jl/blob/master/src/wrap.jl) で検索して、関数の定義を確認する
3. [ArrayFire 本家](http://arrayfire.org/docs/modules.htm) で渡すべき値を確認する

今回使った `convolve2` を例にすると、

```julia
convolve2(signal::AFArray,filter::AFArray,mode::af_conv_mode,domain::af_conv_domain)
```

が関数の定義で、４つの引数（デフォルトも設定されていない）が必要だとわかる。
[convolve2](http://arrayfire.org/docs/group__signal__func__convolve2.htm) を見ると、`signal`, `filter` のほかに、 [convMode](http://arrayfire.org/docs/group__signal__func__convolve2.htm) と [convDomain](http://arrayfire.org/docs/defines_8h.htm#a8519ae1caafab760d2421cd6ac13a3db) を設定すればよいとわかる。
以上をまとめて、次のように呼び出せば GPU での２次元畳込みができる。

```julia
convolve2(s, kernel, AF_CONV_EXPAND, AF_CONV_SPATIAL)
```

## 可視化

状態の経時変化はなんとしても見てみたかった。
既存のライブラリだと Plots.jl がアニメーション機能を有するが、プロットするのに `heatmap` を使用することになる。しかし、可能ならば１セルを１ピクセルに対応させたいという希望があったので、自分でモジュールを書いてみた。

コードは　[Gist](https://gist.github.com/Lirimy/127d6a557ba10d549d17c6687ee9c5c2) においてある。
中身は `animation.jl` をもとにいじっただけだが、速度的にずっと改善してある（３倍くらい？）。
`Plots.jl` では、中間ファイルに PNG を使用しており、実は PNG のエンコーディングに時間がかかっているようなので、中間ファイル形式を BMP に変更した。

[Julia におけるアニメーションのしくみ](https://qiita.com/Lirimy/items/8976bb4dcf5febad178e)

# ArrayFire.jl を用いた実装

## GPU に合わせたコーディング

GPU でやると極端に遅い処理

- 条件分岐・個別の処理
- インデックス操作（読むのはOK、書くのがダメ）

これらをアルゴリズム上うまく回避して書かねばならない。

他の注意としては、特に GeForce などのコンシューマ向け GPU は Float32 に特化しており、うっかり Float64 にキャストされると遅くなってしまう。

```julia
s = rand(AFArray{Float32}, 100, 100)

s / 10   # OK
0.1 * s  # NG

Float32(0.1) * s
0.1f0 * s
```

### 隣接するセルの和

カーネルを畳み込むことで実現した。

```julia
kernel = AFArray(Float32.([0 1 0; 1 0 1; 0 1 0]))
```

### 周期境界条件

行列積で実現した。

[畳み込み時のパディングを工夫してみる](https://qiita.com/Lirimy/items/102c6828f31297bb25ad)

### 比較

比較演算の結果 (Bool) を 1/0 と見て、素直にキャストなしに演算できたのでそれでも良いが、今回は `signbit` を主に利用してみた。

`signbit` は [sign()](http://arrayfire.org/docs/group__arith__func__sign.htm) のラップで、負のときのみ１、他の場合は０を出力する。特筆すべきは、出力が Float32 なので、一貫して同じ型で計算できることである。が、速度にどの程度影響するかは確かめられていない。

```julia
a > b         # Bool
signbit(b-a)  # Float32
```

### 逐次処理

イジングモデルの計算の本来の形は

1. どこか１つのセルに注目して評価する
2. セルの値を更新して、次のセルへ

であるが、GPU だと素直に実装できない。

~~そこで、操作を少々変更して、全体のうちランダムな 1/10 のセルの評価・更新を繰り返すようにした。結果を見るとうまくいっていそうな気はするが、問題があれば教えてほしい。~~

格子模様を描き、色ごとに計算すればよい。
https://qiita.com/Lirimy/items/e25efa016c54d07368c8

## メモリ管理

メモリ管理をせずに計算していると、唐突にこんなエラーが出現する。

```
GPU is out of memory, to avoid this in the future you can:
  @afgc function f(input)   # free all temporary variables inside the function scope
  @afgc a = b + c           # free all temporary arrays inside the assignment scope
  @afgc a .= b + c          # replace with a new value, free the old and all temps
  afgc(threshold)           # garbage collect after GPU memory usage reaches threshold
  finalize(array)           # manually free GPU memory
```

しかし、何も考えず頻繁にガベージコレクションを呼び出しているとパフォーマンスが大きく低下する。
下手なことをやっていたときは　gc (CPU側) が処理時間の 99% を占めていた。
あまりにもひどいので `gc_enable(false)` などと禁断の呪文を呼び出していたが、無事解決策が見つかったので、これからやる人は参考にしてほしい。

### メインループ内で afgc()

手軽さを求めるならこれが良さそう。

```julia
const MB = UInt64(1024*1024);

for i in 1:100 # メインループ
    
    # 処理
    
    afgc(3200MB)
end
```

`nvidia-smi -l 1` でメモリ使用量を見て、できるだけ大きいしきい値を設定しよう。
`afgc()` は CPU/GPU 両方のガベージコレクションを呼び出す。

### 関数定義に　`@afgc` をつける

この方式だと、 `@afgc` が　256 回呼ばれるごとに afgc() が呼ばれる。
上の方式でうまくいったので採用しなかったが、やり方を残しておく。

```julia
@noinline @afgc function update!(s::AFArray{Float32, 2})
    
    # s を更新する
    
    return nothing
end
```

- 戻り値はなし
- 念のため、明示的にスコープをわける

### CPU 側の gc

基本的には `afgc()` 内で呼ばれるので問題はないはずなのだが、画像出力せずにベンチマークしていると、なぜか画像出力したときの数倍遅くなってしまった。よくわからないが、とりあえず `gc(false)` を入れておけば解決はした。false はフルで gc しないという意味らしい。

# まとめ

## 成果

- イジングモデルを大まかに理解できた
- GPU の計算手法を確立できた
    - GPU の様々な制約への対応
    - メモリ管理
- ２次元配列の可視化手法を確立できた

## 課題

- Julia の `gc()` の挙動を理解する
- そもそもゴミを出さない、または減らす手法を確立する
- パッケージの取り扱いを理解する（よく名前が衝突するため）