[iOS] MetalでGPUコンピューティング (2) 群知能


[前回](http://qiita.com/yuky_az/items/ece9b64befc635e89f1a)はMetalによるGPUコンピューティング(GPGPU)の基礎と特性を扱いましたが、今回はその応用を扱います。アプリの描画部分にMetalを用いるのではなく、ロジック部分にMetalを用いるという通常と逆の試みを行います。
今回は、Metalを用いて[群知能](https://ja.wikipedia.org/wiki/%E7%BE%A4%E7%9F%A5%E8%83%BD)を実装します。
>各個体はローカルに互いと、そして彼らの環境と対話する。個々のエージェントがどう行動すべきかを命じている集中的な制御構造は通常存在しないが、そのようなエージェント間の局所相互作用はしばしば全体の行動の創発(emergence)をもたらす。このようなシステムの自然界の例として、アリの巣、鳥の群れ、動物の群れ、細菌のコロニー、魚の群れなどがある。
>-群知能- Wikipedia

コンピュータ上の群知能として有名なものに、クレイグ・レイノルズによって考案された[Boids](https://ja.wikipedia.org/wiki/%E3%83%9C%E3%82%A4%E3%83%89_(%E4%BA%BA%E5%B7%A5%E7%94%9F%E5%91%BD))があります。
Boidsは複雑系や自己組織化に深く関わっているのですが、各個体が以下のたった3つのルールに従うだけで、個体の集合が魚や鳥の群れのように振る舞います。
1.分離（Separation）各個体は他の個体とぶつからないように動く。
2.整列（Alignment）各個体は他の個体と速度と方向を合わせるように動く。
3.結合（Cohesion）各個体は他の個体の集合の中心方向へ向かうように動く。

上記と全く同じにしても面白くないので、今回はこのルールに少しアレンジを加えました。
a.等距離 各個体が等距離を保つように動く。
b.並行 各個体が他の個体と同じ向きになるように動く。
c.等速 各個体は他の個体と同じ速度になるように動く。

Boidsとの違いは群れが中心を持たず、各個体がより分散的に動くことです。

それでは、上記のa、b、cを数式化してみます。
まず、a.等距離の式からです。
他の個体が近い場合は他の個体から遠ざかる方向に、遠い場合は近ずく方向に角度が変化します。

$ \Delta\Phi_d = \alpha \left(\sum_{k=0}^n (\theta_k - \rho) exp(-b (d_k - a)^2)  - \sum_{k=0}^n (\psi_k - \rho) exp(-b d_k^2) \right) \Delta t$

$\Delta\Phi_d$: 角度の変化
$\alpha$: 角度変化の大きさを決める定数
$n$: 自身を除いた個体数
$\theta_k$: 自身から他の個体へ向かうベクトルの向き(角度)
$\psi_k$: 他の個体から自身へ向かうベクトルの向き(角度)
$\rho$: 自身の向き(角度)
$a$: 個体間の望ましい距離を決める定数
$b$: 影響範囲の広さを決める定数
$d_k$: 自身と他の個体間の距離
$\Delta t$: 時間間隔

右辺()内のの第一項は引力のように個体同士を近づける項、第二項は斥力のように個体同士を遠ざける項です。各項は個体間の距離に依存しており、$a/2$の距離で角度変化は釣り合って$0$になります。

次に、b.並行の式からです。
他の個体と平行になるように、向きを変更します。近い個体ほど、与える影響が大きくなります。

$ \Delta\Phi_p = \beta \sum_{k=0}^n (\omega_k - \rho) exp(-b d_k^2) \Delta t$

$\Delta\Phi_p$: 角度の変化
$\beta$: 角度変化の大きさを決める定数
$n$: 自身を除いた個体数
$\omega_k$: 他の個体の向き(角度)
$\rho$: 自身の向き(角度)
$a$: 個体間の距離を決める定数
$b$: 影響範囲の広さを決める定数
$d_k$: 自身と他の個体間の距離
$\Delta t$: 時間間隔

そして、c.等速の式です。
近い個体との速度差を解消するように、速度が変化します。

$ \Delta V = \gamma \sum_{k=0}^n (v_k-v) exp(-b d_k^2) \Delta t$

$\Delta V$: 速度の変化
$\gamma$: 速度変化の大きさを決める定数
$n$: 自身を除いた個体数
$v_k$: 他の個体の速度
$v$: 自身の速度
$a$: 個体間の距離を決める定数
$b$: 影響範囲の広さを決める定数
$d_k$: 自身と他の個体間の距離
$\Delta t$: 時間間隔

それでは、コードの方を記述します。
ロジック部分にはMetalを、描画部分はUIImageViewを使いました。
本来はSpriteKitを使用したかったのですが、SpriteKitを用いるとバックで動作するMetalが干渉してしまうのか、画面が頻繁にちらつく問題が発生してしまいました。

まず、ViewControllerからです。

```Swift
import UIKit
import Metal
import SpriteKit

class ViewController: UIViewController {

    let nodeCount = 200
    var nodes:[Node] = []
    
    var metalController:MetalController!
    
    var displayLink:CADisplayLink!
    var lastTimeStamp:CFTimeInterval = 0
    
    var fishes:[Fish] = []
    
    override func viewDidLoad()
    {
        super.viewDidLoad()
        
        let width = self.view.frame.size.width;
        let height = self.view.frame.size.height;
        
        for _ in 0..<nodeCount {
            let randX = Float(arc4random_uniform(UInt32(width)))
            let randY = Float(arc4random_uniform(UInt32(height)))
            let randVX = Float(arc4random_uniform(UInt32(100)))/100.0
            let randVY = Float(arc4random_uniform(UInt32(100)))/100.0
            let randAngle = Float(arc4random_uniform(UInt32(100)))/100.0 * Float(M_PI) * 2 - Float(M_PI)

            nodes.append(Node(positionX: randX, positionY: randY, velocityX: randVX, velocityY: randVY, angle: randAngle))
            let fish = Fish()
            fish.center = CGPoint(x: CGFloat(randX), y: height-CGFloat(randY))
            self.view.addSubview(fish)
            fishes.append(fish)
        }

        metalController = MetalController(nodes, width: width, height: height)
        
        displayLink = CADisplayLink(target: self, selector: #selector(ViewController.update))
        displayLink.add(to: RunLoop.current, forMode: .defaultRunLoopMode)
    }
    
    func update()
    {
        if lastTimeStamp == 0 {
            lastTimeStamp = displayLink.timestamp
            return
        }
        
        let now = displayLink.timestamp
        let interval = now - lastTimeStamp
        lastTimeStamp = now
        
        metalController.move(nodes: nodes, interval: Float(interval)){( resultNodes:[Node]) -> Void in
            nodes = resultNodes
            for (i, fish) in fishes.enumerated() {
                fish.setNode(node: nodes[i], height: self.view.frame.size.height)
            }
        }
    }
    
    deinit {
        displayLink.invalidate()
    }
}

struct Node
{
    var positionX: Float = 0
    var positionY: Float = 0
    var velocityX: Float = 0
    var velocityY: Float = 0
    var angle: Float = 0
}
```
次に、UImageViewを継承した描画用のクラスです。

```Swift
import UIKit

class Fish: UIImageView {

    init(){
        super.init(frame: CGRect(x: 0, y: 0, width: 30, height: 30))
        self.image = UIImage(named: "fish.png")
    }
    
    func  setNode(node:Node, height:CGFloat) {
        self.center = CGPoint(x: CGFloat(node.positionX), y: height - CGFloat(node.positionY))
        self.transform =  CGAffineTransform(rotationAngle: CGFloat(node.angle))
    }
    
    required init?(coder aDecoder: NSCoder) {
        fatalError("init(coder:) has not been implemented")
    }
}
```
そして、Metalによる演算を行うためのクラスです。

```Swift
import UIKit
import Metal

class MetalController: NSObject {
    
    private var device: MTLDevice!
    private var defaultLibrary: MTLLibrary!
    private var commandQueue: MTLCommandQueue!
    private var computePipelineState: MTLComputePipelineState!
    
    private var threadsPerThreadgroup:MTLSize!
    private var threadgroupsCount:MTLSize!
    
    private var nodeCountBuffer: MTLBuffer!
    private var widthBuffer: MTLBuffer!
    private var heightBuffer: MTLBuffer!
    private var outBuffer: MTLBuffer!

    let nodeCount:UInt32
        
    init(_ nodes:[Node], width:CGFloat, height:CGFloat)
    {
        var count = UInt32(nodes.count)
        self.nodeCount = count
        
        device = MTLCreateSystemDefaultDevice()
        defaultLibrary = device.newDefaultLibrary()
        commandQueue = device.makeCommandQueue()
        let ml2Func = defaultLibrary.makeFunction(name: "move")!
        computePipelineState = try! device.makeComputePipelineState(function: ml2Func)
        
        let threadWidth = 64
        threadsPerThreadgroup = MTLSize(width: threadWidth, height: 1, depth: 1)
        threadgroupsCount = MTLSize(width: (Int(nodeCount) + threadWidth - 1) / threadWidth, height: 1, depth: 1)
        
        nodeCountBuffer = device.makeBuffer(bytes: &count, length: MemoryLayout.size(ofValue: count), options: [])
        var wdth = Float(width)
        widthBuffer = device.makeBuffer(bytes: &wdth, length: MemoryLayout.size(ofValue: wdth), options: [])
        var hght = Float(height)
        heightBuffer = device.makeBuffer(bytes: &hght, length: MemoryLayout.size(ofValue: hght), options: [])
        outBuffer = device.makeBuffer(bytes: nodes, length: nodes.byteLength, options: [])
    }
    
    func move( nodes:[Node], interval:Float, callBack:([Node]) -> Void)
    {
        let commandBuffer = commandQueue.makeCommandBuffer()
        let computeCommandEncoder = commandBuffer.makeComputeCommandEncoder()
        computeCommandEncoder.setComputePipelineState(computePipelineState)
        
        let inBuffer = device.makeBuffer(bytes: nodes, length: nodes.byteLength, options: [])
        var itvl = Float32(interval)
        let intervalBuffer = device.makeBuffer(bytes: &itvl, length: MemoryLayout.size(ofValue: itvl), options: [])
        
        computeCommandEncoder.setBuffer(inBuffer, offset: 0, at: 0)
        computeCommandEncoder.setBuffer(nodeCountBuffer, offset: 0, at: 1)
        computeCommandEncoder.setBuffer(intervalBuffer, offset: 0, at: 2)
        computeCommandEncoder.setBuffer(widthBuffer, offset: 0, at: 3)
        computeCommandEncoder.setBuffer(heightBuffer, offset: 0, at: 4)
        computeCommandEncoder.setBuffer(outBuffer, offset: 0, at: 5)
        
        computeCommandEncoder.dispatchThreadgroups(threadgroupsCount, threadsPerThreadgroup: threadsPerThreadgroup)
        computeCommandEncoder.endEncoding()
        
        commandBuffer.commit()
        commandBuffer.waitUntilCompleted()
        
        let data = Data(bytesNoCopy: outBuffer.contents(), count: nodes.byteLength, deallocator: .none)
        var resultNodes = [Node](repeating: Node(positionX: 0, positionY: 0, velocityX: 0, velocityY: 0, angle: 0), count: nodes.count)
        resultNodes = data.withUnsafeBytes {
            Array(UnsafeBufferPointer<Node>(start: $0, count: data.count/MemoryLayout<Node>.size))
        }

        callBack(resultNodes)
    }

}

private extension Array {
    var byteLength: Int {
        return self.count * MemoryLayout.size(ofValue: self[0])
    }
}

```
最後に、GPU側のコード、シェーダーです。上記の数式は、この中でコードにします。

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

constant float alpha = 0.025;
constant float beta = 0.1;
constant float gamma = 0.005;

constant float spaceRatio = 0.12;

struct Node
{
    float positionX;
    float positionY;
    float velocityX;
    float velocityY;
    float angle;
};

static float getDistace(float x1, float y1, float x2, float y2)
{
    float dx = x1-x2;
    float dy = y1-y2;
    return sqrt(dx*dx + dy*dy);
}

//From -pi to pi
static float getRangedAngle(float angle)
{
    if (angle > M_PI_F){
        angle -= 2 * M_PI_F;
    }else if (angle < -M_PI_F){
        angle += 2 * M_PI_F;
    }
    return angle;
}

kernel void move(const device Node *inNode [[ buffer(0) ]],
                 const device uint &nodeCount [[ buffer(1) ]],
                 const device float &interval [[ buffer(2) ]],
                 const device float &width [[ buffer(3) ]],
                 const device float &height [[ buffer(4) ]],
                 device Node *outNode [[ buffer(5) ]],
                 uint id [[ thread_position_in_grid ]])
{
    Node currentNode = inNode[id];
    
    float a = width * spaceRatio;
    float b = 6.25 / a / a;
    float dAngle = 0;
    
    float velocityX = currentNode.velocityX;
    float velocityY = currentNode.velocityY;
    float velocity = sqrt(velocityX*velocityX + velocityY*velocityY);
    
    float outerSpace = width * 0.1;
    
    for (uint i=0; i<nodeCount; i++){
        if (i == id){
            continue;
        };
        
        Node node = inNode[i];
        
        float distance = getDistace(node.positionX, node.positionY, currentNode.positionX, currentNode.positionY);
        
        float nearAngle = getRangedAngle(atan2(node.positionX-currentNode.positionX, node.positionY-currentNode.positionY) - currentNode.angle);
        float farAngle = getRangedAngle(atan2(currentNode.positionX-node.positionX, currentNode.positionY-node.positionY) - currentNode.angle);
        float attraction = exp(-b * (distance - a)*(distance - a));
        float repulsion = exp(-b * distance * distance);
        dAngle += alpha * (nearAngle*attraction + farAngle*repulsion)*interval;
        
        float parallelAngleDif = getRangedAngle(node.angle - currentNode.angle);
        dAngle += beta * parallelAngleDif * exp(-b * distance * distance) * interval;
        
        float nodeVelocity = sqrt(node.velocityX*node.velocityX + node.velocityY*node.velocityY);
        velocity += gamma * (nodeVelocity - velocity) * exp(-b * distance * distance);;
    }
    
    float newAngle = getRangedAngle(currentNode.angle + dAngle);
    velocityX = velocity * sin(newAngle);
    velocityY = velocity * cos(newAngle);
    
    float newPositionX = currentNode.positionX + velocityX;
    if (newPositionX > (width + outerSpace)) {
        newPositionX -= (width + outerSpace*2);
    }
    if (newPositionX < -outerSpace) {
        newPositionX += width + outerSpace*2;
    }
    
    float newPositionY = currentNode.positionY + velocityY;
    if (newPositionY > height + outerSpace) {
        newPositionY -= height + outerSpace*2;
    }
    if (newPositionY < -outerSpace) {
        newPositionY += height + outerSpace*2;
    }
    
    outNode[id].positionX = newPositionX;
    outNode[id].positionY = newPositionY;
    outNode[id].velocityX = velocityX;
    outNode[id].velocityY = velocityY;
    outNode[id].angle = newAngle;
}
```
それでは、動作を見てみます。検証機にはiPhone6 Plusを用いました。
以下は$\alpha = 0.03$、$\beta = -0.35$、$\gamma = 0.005$の場合です。
下の画像をクリックすると動画が再生します。
 [![群知能1](http://img.youtube.com/vi/cAUOJ-wmeUA/0.jpg)](http://www.youtube.com/watch?v=cAUOJ-wmeUA)
最初はそれぞれの個体がランダムな方向を向いており、ランダムな速度です。
![スクリーンショット 2016-11-10 22.35.36.png](https://qiita-image-store.s3.amazonaws.com/0/35297/dd1df74f-0d3e-f3e2-a513-daa2f65d7068.png "スクリーンショット 2016-11-10 22.35.36.png")
時間が経過すると、多くの魚が群れを形成し同じ方向に動く傾向が見えてきます。また、小規模なグループが一斉に旋回する様子も観察できます。イメージとしては魚の群れや鳥の群れに近いかと思います。
![スクリーンショット 2016-11-10 22.59.46.png](https://qiita-image-store.s3.amazonaws.com/0/35297/552d2254-1584-4b66-f5c4-5e293e225f89.png "スクリーンショット 2016-11-10 22.59.46.png")


以下は$\alpha = 0.16$、$\beta = 0$、$\gamma = 0.005$の場合です。
下の画像をクリックすると動画が再生します。

[![群知能2](http://img.youtube.com/vi/GjdrYERYKVc/0.jpg)](http://www.youtube.com/watch?v=GjdrYERYKVc)
この場合、魚たちが輪になって回転を始めます。内外二つの輪がそれぞれ逆方向に回転します。歯車や自動車のホイールのようです。
![スクリーンショット 2016-11-10 23.23.05.png](https://qiita-image-store.s3.amazonaws.com/0/35297/197685e4-e1c0-8587-cc57-34045df831de.png "スクリーンショット 2016-11-10 23.23.05.png")

そして、これが今回最も面白かったのですが、以下は$\alpha = 0.025$、$\beta = 0.2$、$\gamma = 0.005$の場合です。
下の画像をクリックすると動画が再生します。

[![群知能3](http://img.youtube.com/vi/WlMiJDzC2f4/0.jpg)](http://www.youtube.com/watch?v=WlMiJDzC2f4)
魚たちが集まって、まるで一つの生き物であるような動きをします。姿形はクラゲのようでありますが、無数の個体が集まって一つの個体を形成する様子は粘菌のようでもあり、魚を細胞とするならば多細胞生物に例えることもできるかと思います。
![スクリーンショット 2016-11-10 23.45.36.png](https://qiita-image-store.s3.amazonaws.com/0/35297/a515c201-3d26-8f88-f21e-cfcf5151c861.png "スクリーンショット 2016-11-10 23.45.36.png")

今回は、Metalを純粋にロジックのみに用いてアプリを実装してみました。
200の個体を用いたのですが、それぞれの相互作用を考えたので200x200で40000回のの演算を毎フレームごと、1秒間に60回行ったことになります。しかも、三角関数、指数、平方根などを含むCPUではコストの高い処理を多く含んでいました。あまりパフォーマンスを考慮したコードになっていないので、無駄な処理も多く含まれています。
それにも関わらず、フレームが落ちることは全くありませんでした。なおかつ、CPUの消費はわずか35%ほどでした。
exp関数だけでも秒間1000万回（！）使用していることになります。

改めて、MetalによるGPUコンピューティングの威力を感じた次第です。今回は群知能を対象としましたが、もちろん機械学習などへの応用も期待できるかと思います。

今回使用したコードは、こちらに置いておきます。
[https://github.com/yukinaga/SwarmIntelligence](https://github.com/yukinaga/SwarmIntelligence)
パラメータを変更するともっと面白い世界が見られるかもしれないので、ぜひ試してくださいね。