Metalシェーダによるマルチエージェントシミュレーション

### Metalとは
MetalはiOS8で登場したローレベルなグラフィックAPIで、画像処理や並列演算を得意とします。
このMetalによる並列演算を利用してマルチエージェントシミュレーションを試してみたいと思います。

### マルチエージェントシステムとは
マルチエージェントシミュレーションとは、複数のエージェント（人や生物など）で集団を形成し、行動ルールのもとに相互作用をシミュレーションするシステムです。

例えば以下の図のように個人（エージェント）を見ると「進む」というルールに従って行動していますが、集団でみると、近くの人に「追従」するというルールで行動しています。
この相互作用が積み重なると、みんなが同じ方向に進むという結果になります。
<img src="https://qiita-image-store.s3.amazonaws.com/0/14751/2a8cfa90-882c-c3a7-a30e-cd8a6c9fe12a.png" width=400>

このような相互作用は「鳥の群れ」「道路の渋滞」「流行」など身近にも存在しています。
マルチエージェントシミュレーションは、これらの現象についてその仕組みを解析するのに用いられます。

今回は "Dynamic Models of Segregation"（住み分けの動学的モデル）のシミュレーションを行いました。

たとえ白人と黒人が隣同士で暮らすことに抵抗がなくとも、いつの間にか白人が多く居住する地域と黒人が多く居住する地域に分かれてしまう理由を考察するためのモデルで、1971年に[Thomas Crombie Schelling](https://ja.wikipedia.org/wiki/%E3%83%88%E3%83%BC%E3%83%9E%E3%82%B9%E3%83%BB%E3%82%B7%E3%82%A7%E3%83%AA%E3%83%B3%E3%82%B0)によって発表されました。

ルールは以下の通りです。

> 2種類の人種が存在していて最初はランダムに配置されます。
> 各人は自分の周囲にいる人のうち同じ人種が50%を満たしていれば満足してその場にとどまり、50%に満たない場合はその場を移動します。


これを繰り返すことでどのような結果が得られるかをシミュレーションします。

### 実装
それではやっとコードの話に進みます。
先にサンプルプロジェクトのリンクを貼っておきます。[サンプル](https://github.com/yoshimin/Segregation)
ところどころ省略されているので分かりづらかったらサンプルを動かしてみてください。

今回のマルチエージェントシミュレーションを実行するための流れは以下のようになります。

1. 各エージェントをオブジェクトと `UIView`をエージェントの数分用意。
2. 生成したオブジェクトの配列をシェーダに渡し、ルールに従って座標を計算。
3. `UIView` に計算した座標を反映。
4. 以降繰り返し。

#### 準備
Metalで並列演算を動かすための初期化処理は5ステップです。
以下のものを作成しています。

1. エージェントを表すオブジェクト
2. MTLDevice
3. Metalシェーダ
4. MTLComputePipelineState
5. MTLCommandQueue

##### 1. エージェントを表すオブジェクト

以下のようなstructで表現することとします。
エージェントの数だけ生成しておきます。

```
struct Object {
    var group: UInt32 = 0
    var x: Float = 0
    var y: Float = 0
    var angle: Float = 0
}
```
angleはエージェントが移動する際の方向として使用します。

##### 2. MTLDevice
まずは当たり前ですがMetalフレームワークをimportします。

```
import Metal
```

`MTLDevice` は `MTLCreateSystemDefaultDevice` によって取得できます。

```
var device: MTLDevice = MTLCreateSystemDefaultDevice()!
```

##### 3. Metalシェーダ
Metalシェーダは Metal Shading Language という、C++をベースとする独自言語を使用します。

```
#include <metal_stdlib>
using namespace metal;

struct Object {
    uint group;
    float x;
    float y;
    float angle;
};

float getDistance(float2 pos1, float2 pos2) {
    float2 d = pos1 - pos2;
    return sqrt(d.x*d.x+d.y*d.y);
}

kernel void simulate(const device Object *in [[ buffer(0) ]],
                     const device uint &count [[ buffer(1) ]],
                     const device float &width [[ buffer(2) ]],
                     const device float &height [[ buffer(3) ]],
                     device Object *out [[ buffer(4) ]],
                     uint id [[ thread_position_in_grid ]]) {
    Object object = in[id];
    out[id] = object;

    float2 objectPos = float2(object.x,object.y);

    int total = 0;
    int group = 0;
    for (uint i=0; i<count; i++){
        if (i == id){
            continue;
        };

        Object other = in[i];

        float2 otherPos = float2(other.x,other.y);
        float dist = getDistance(objectPos,otherPos);

        int near = step(dist,20);
        total += near;
        group += (object.group == other.group) * near;
    }
    
    int stay = step(1.0,total)*step(total*0.5, group);
    if (stay) {
        return;
    }
    
    float translationX = sin(object.angle);
    float translationY = cos(object.angle);

    float x = object.x + translationX;
    x -= (x > width) ? width : 0;
    x += (x < 0) ? width : 0;

    float y = object.y + translationY;
    y -= (y > height) ? height : 0;
    y += (y < 0) ? height : 0;

    out[id].x = x;
    out[id].y = y;
}
```
`simulate` 関数の引数は以下のとおりです。
第一引数: `Object in` エージェントの配列
第二引数: `uint count` エージェントの数
第二引数: `float width` 画面の幅
第二引数: `float height` 画面の高さ
第二引数: `Object out` 結果の出力先
第二引数: `uint id` エージェントのインデックス
 `thread_position_in_grid` でエージェントのインデックスが取得できるので、各エージェントは `in[id]` で取得できます。
_※ `thread_position_in_grid` についてはまたあとで説明が出てきます。_
あとは先述のルールに従って計算処理を書いていき、 `out` に計算結果を反映させます。

##### 4. MTLComputePipelineState

Metalシェーダを作成したので、これをコンパイルする必要があります。
`MTLComputePipelineState` プロトコルは、コンパイルされたプログラムへの参照をエンコードするために使用される軽量オブジェクトのインターフェイスを定義します。 
`MTLComputePipelineState` の生成は負荷が高いのでメンバ変数として定義し初期化時に生成します。

```
var pipelineState: MTLComputePipelineState!
```

```
let defaultLibrary = device.makeDefaultLibrary()
let function = defaultLibrary.makeFunction(name: "simulate")!
pipelineState = try! device.makeComputePipelineState(function: function)
```

`device.makeDefaultLibrary()` で取得した `MTLLibrary` オブジェクトを使ってコンパイルしたシェーダへアクセスすることができます。
シェーダの関数名を指定し `MTLFunction` を作成したら、 `device.makeComputePipelineState(function:)` で `MTLFunction` をコンパイルしたコードに変換します。


##### 5. MTLCommandQueue

最後に `MTLCommandQueue` を生成します。
`MTLCommandQueue` プロトコルは、GPUが実行するコマンド・バッファの順序付きリストをキューに入れることができるオブジェクトのインタフェースを定義します。

これもメンバ変数として定義しておきます。

```
var commandQueue: MTLCommandQueue!
```

```
commandQueue = device.makeCommandQueue()
```

#### 演算処理の実行
並列演算を動かすための準備ができたので次は実行するための処理を書いていきます。
実行用に `func execute() -> [Object]` というメソッドを用意して以下の処理を書いていきます。

1. `MTLCommandBuffer` の作成
2. 実行コマンドをエンコード
3. `MTLCommandBuffer` をコミット
4. 計算結果を返す

##### 1. MTLCommandBuffer の作成
実行したいコマンドのリストとなります。

```
let commandBuffer = commandQueue?.makeCommandBuffer()
```

##### 2. 実行コマンドをエンコード

```
let commandBuffer = commandQueue.makeCommandBuffer()

let encoder = commandBuffer?.makeComputeCommandEncoder()
encoder?.setComputePipelineState(pipelineState)

let buffer = device?.makeBuffer(bytes: objects, length: objects.byteLength, options: [])
var count = UInt32(objects.count)
let countBuffer = device?.makeBuffer(bytes: &count, length: MemoryLayout.size(ofValue: count), options: [])
var width = Float(UIScreen.main.bounds.size.width)
var height = Float(UIScreen.main.bounds.size.height)
let widthBuffer = device?.makeBuffer(bytes: &width, length: MemoryLayout.size(ofValue: width), options: [])
letheightBuffer = device?.makeBuffer(bytes: &height, length: MemoryLayout.size(ofValue: height), options: [])
let outBuffer = device?.makeBuffer(bytes: objects, length: objects.byteLength, options: [])!

encoder?.setBuffer(buffer, offset: 0, index: 0)
encoder?.setBuffer(countBuffer, offset: 0, index: 1)
encoder?.setBuffer(widthBuffer, offset: 0, index: 2)
encoder?.setBuffer(heightBuffer, offset: 0, index: 3)
encoder?.setBuffer(outBuffer, offset: 0, index: 4)

let groupsize = MTLSize(width: 64, height: 1, depth: 1)
let numgroups = MTLSize(width: (objects.count + groupsize.width - 1) / groupsize.width, height: 1, depth: 1)

encoder?.dispatchThreadgroups(numgroups, threadsPerThreadgroup: groupsize)
encoder?.endEncoding()
```

コマンドをエンコードするためのエンコーダ（`MTLComputeCommandEncoder`）を作成します。
`MTLComputeCommandEncoder` オブジェクトを作成したら、このオブジェクトを使用して、並列演算処理コマンドをエンコードします。
先程作成した、実行される関数を含む `MTLComputePipelineState` オブジェクトをセットします。
入力データと出力先を保持するリソースを指定します。（指定するリソースは `simulate` 関数の引数で定義した5つです。）
`dispatchThreadgroups(_: threadsPerThreadgroup:)` メソッドを呼び出して、グリッド用の指定された数のスレッドグループとスレッドグループあたりのスレッド数で計算関数をエンコードします。

Metalシェーダの項で `thread_position_in_grid` が登場しましたが、`dispatchThreadgroups(_: threadsPerThreadgroup:)` で指定するサイズによってグリッドの範囲が決まり、`thread_position_in_grid` でそのグリッド内のインデックスが取得できるということだそうです。

`endEncoding()`を呼び出してエンコードを終了します。

##### 3. MTLCommandBuffer をコミット
コマンドを実行します。

```
commandBuffer?.commit()
commandBuffer?.waitUntilCompleted()
```

#### 4. 計算結果を返す
3で実行すると演算処理が走り、outBufferから結果を取得し返します。

```
let data = Data(bytesNoCopy: outBuffer.contents(), count: objects.byteLength, deallocator: .none)
objects = data.withUnsafeBytes {
    [Object](UnsafeBufferPointer<Object>(start: $0, count: data.count/MemoryLayout<Object>.size))
}
return objects
```

あとは `UIView` 計算した座標を反映させれば完了です。

### 実行結果
最初はランダムに配置されていますが、シミュレートを重ねると右のように同じ人種同士で固まるという結果が得られました。

<img width="200" alt="IMG_7120.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/14751/f3433a93-c1c1-6577-364e-1140b8939f00.png"> → <img width="200" alt="IMG_7117.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/14751/fe1f1d1d-7ca4-0ba0-655d-d4cd62e69ff7.png">

周囲の50%を同人種で満たしていればOKということは残り約半数は異人種でもいいということになるので結構寛容的に思われますが、それでもこのようにはっきり同人種同士で固まることになるのですね。
試しに許容の水準を30%、10%に下げてみてもある程度分居される結果になりました。

|30%|10%|
|---|---|
|<img width="120" alt="IMG_7122.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/14751/860f259e-8ee8-8812-ce51-04b96db83632.png">|<img width="120" alt="IMG_7124.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/14751/db0e8f58-dafa-b93a-fcd3-02c508b4af00.png">
|

シェリングは現象を隣人に対する寛容性あるいは我慢の強弱に起因するものであると仮定していますが、確かに個々が寛容的であっても分居は発生し、そのレベルは許容の水準によるものが影響しているようにみえます。

### ネイティブでやった場合の比較
最後に、Metalを使うとどれくらい良いのか、ネイティブとの比較を載せておきます。

#### 実行速度
60回シミュレーションを実行した平均はそれぞれ以下の通りです

|エージェント数|Metal|ネイティブ|
|---|---|---|
|100|0.002582641443|0.007742651304|
|500|0.003791083892|0.03941152493|
|1000|0.004310701291|0.1502834062|
|1500|0.005150139332|0.3368620038|
|2000|0.00574699839|0.5972672562|

#### CPU
エージェント数2000のCPUです

|Metal|ネイティブ|
|---|---|
|<img width="400" alt="スクリーンショット 2018-12-10 23.51.59.png" src="https://qiita-image-store.s3.amazonaws.com/0/14751/ec867d69-f9d1-0478-1426-b0dc3f627d65.png">|<img width="400" alt="スクリーンショット 2018-12-10 23.53.21.png" src="https://qiita-image-store.s3.amazonaws.com/0/14751/06b124a1-4b63-e26c-c0c4-a4ef4a59ede9.png">|

どちらも圧倒的にMetalのほうがパフォーマンスが良いですね。

### まとめ
マルチエージェントシミュレーションをMetalで試してみましたが、Metalシェーダ（Metal Shading Language）を書くというハードルはありますが、Metalフレームワークを使うと割と簡単に動かせました。
マルチエージェントシミュレーションについては今回は簡単なモデルで試しましたが、他にもいろいろモデルがるのでシェーダの練習がてらシミュレーションしてみるのも良いかもしれません。
自分でも身近にある現象と理由を考察してシミュレーションしてみると面白そうですね。（私にはハードル高いけど...）

### 参考
https://qiita.com/yuky_az/items/3f3de1034e0600be490d
http://mas.kke.co.jp/modules/tinyd4/index.php?id=13