Xcode
iOS
Metal
GPU
Swift

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

More than 1 year has passed since last update.

前回はMetalによるGPUコンピューティング(GPGPU)の基礎と特性を扱いましたが、今回はその応用を扱います。アプリの描画部分にMetalを用いるのではなく、ロジック部分にMetalを用いるという通常と逆の試みを行います。

今回は、Metalを用いて群知能を実装します。


各個体はローカルに互いと、そして彼らの環境と対話する。個々のエージェントがどう行動すべきかを命じている集中的な制御構造は通常存在しないが、そのようなエージェント間の局所相互作用はしばしば全体の行動の創発(emergence)をもたらす。このようなシステムの自然界の例として、アリの巣、鳥の群れ、動物の群れ、細菌のコロニー、魚の群れなどがある。

-群知能- Wikipedia


コンピュータ上の群知能として有名なものに、クレイグ・レイノルズによって考案されたBoidsがあります。

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からです。

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を継承した描画用のクラスです。

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による演算を行うためのクラスです。

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側のコード、シェーダーです。上記の数式は、この中でコードにします。

#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

最初はそれぞれの個体がランダムな方向を向いており、ランダムな速度です。

スクリーンショット 2016-11-10 22.35.36.png

時間が経過すると、多くの魚が群れを形成し同じ方向に動く傾向が見えてきます。また、小規模なグループが一斉に旋回する様子も観察できます。イメージとしては魚の群れや鳥の群れに近いかと思います。

スクリーンショット 2016-11-10 22.59.46.png

以下は$\alpha = 0.16$、$\beta = 0$、$\gamma = 0.005$の場合です。

下の画像をクリックすると動画が再生します。

群知能2

この場合、魚たちが輪になって回転を始めます。内外二つの輪がそれぞれ逆方向に回転します。歯車や自動車のホイールのようです。

スクリーンショット 2016-11-10 23.23.05.png

そして、これが今回最も面白かったのですが、以下は$\alpha = 0.025$、$\beta = 0.2$、$\gamma = 0.005$の場合です。

下の画像をクリックすると動画が再生します。

群知能3

魚たちが集まって、まるで一つの生き物であるような動きをします。姿形はクラゲのようでありますが、無数の個体が集まって一つの個体を形成する様子は粘菌のようでもあり、魚を細胞とするならば多細胞生物に例えることもできるかと思います。

スクリーンショット 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

パラメータを変更するともっと面白い世界が見られるかもしれないので、ぜひ試してくださいね。