【ALife】OpenGLで爆速シミュレーション

はじめに

　OpenGL4.0が発表されてから10年近く経ちました．その間にOpenGLよりも低レベルなAPIであるVulkanが発表され，3DグラフィックスAPIの低レベル化が決定的なものとなったわけですが，グラフィックスリソースの初期化や管理などの煩雑な処理を自動化してくれる上位のAPIとしてOpenGLは一定の価値を残しています．そんなOpenGLの現代的な機能としてプログラマブルシェーダが挙げられると思います．
　本記事では人工生命分野のモデルを題材として，プログラマブルシェーダの中でもあまり目立つことのないジオメトリシェーダとコンピュートシェーダに着目した，CPUのみでは実現できない速度で動作するセルラーオートマトンやマルチエージェントモデルを実装し紹介したいと思います．なお本記事で扱うモデルは「作って動かすALife　実装を通した人工生命モデル理論入門（岡 瑞起，池上 高志，ドミニク・チェン，青木 竜太，丸山 典宏）」を参考にしています．人工生命分野に興味のある方はぜひご覧ください．

　ちなみに今回扱うコードのうち，人工生命モデルのコードはすべてこちらに上がっています．

OpenGLとは

　この項では今回使う技術も合わせてOpenGLの概要に触れていきます．基礎的なことなので飛ばしても問題ありません．

　OpenGLはクロノス・グループが策定しているクロスプラットフォームな3DグラフィックスAPIで，GPUを使うため最適化すれば非常に高速に複雑なシーンを描画することが可能です．また，C/C++だけでなく様々な言語でのバインディングが存在するため，自分が慣れ親しんでいる言語でOpenGLを利用したアプリケーションを開発することができます．今回は簡単にライブラリの導入ができ，ソースコードも読みやすいPythonを使用します．

OpenGLで三角形描画

　現代的なOpenGLでのプログラミングに慣れていただくために，簡単な例として三角形描画を紹介します．

triangle.py

from ctypes import sizeof
import sys

import numpy as np

from OpenGL.GL import *
import glfw

class Shader:
    def __init__(self):
        self.handle = glCreateProgram()
    
    def attach_shader(self, content, type):
        shader = glCreateShader(type)
        glShaderSource(shader, [content])
        glCompileShader(shader)

        status = ctypes.c_uint(GL_UNSIGNED_INT)
        glGetShaderiv(shader, GL_COMPILE_STATUS, status)
        if not status:
            print(glGetShaderInfoLog(shader).decode("utf-8"), file=sys.stderr)
            glDeleteShader(shader)
            return False
        
        glAttachShader(self.handle, shader)
        glDeleteShader(shader)
        return True

    def link(self):
        glLinkProgram(self.handle)
        status = ctypes.c_uint(GL_UNSIGNED_INT)
        glGetProgramiv(self.handle, GL_LINK_STATUS, status)
        if not status:
            print(glGetProgramInfoLog(self.handle).decode("utf-8"), file=sys.stderr)
            return False
        return True
    
    def use(self):
        glUseProgram(self.handle)

    def unuse(self):
        glUseProgram(0)

vert = """
#version 450

layout (location = 0) in vec2 vertex;
layout (location = 1) in vec3 color;

layout (location = 0) out vec3 outColor;

void main() {
    gl_Position = vec4(vertex, 0, 1);
    outColor = color;
}
"""

frag = """
#version 450

layout (location = 0) in vec3 color;

layout (location = 0) out vec4 outColor;

void main() {
    outColor = vec4(color, 1);
}
"""

def main():
    if not glfw.init():
        raise RuntimeError("failed to initialize GLFW")

    window = glfw.create_window(700, 700, "triangle", None, None)
    if not window:
        glfw.terminate()
        raise RuntimeError("failed to create GLFWwindow")
    
    glfw.window_hint(glfw.CONTEXT_VERSION_MAJOR, 4)
    glfw.window_hint(glfw.CONTEXT_VERSION_MINOR, 5)
    glfw.make_context_current(window)

    # シェーダの生成
    program = Shader()
    program.attach_shader(vert, GL_VERTEX_SHADER)
    program.attach_shader(frag, GL_FRAGMENT_SHADER)
    program.link()

    data = np.array([
        0, 1,  1, 0, 0, # 頂点(0, 1)，  色(1, 0, 0)
        1, -1,  0, 1, 0,# 頂点(1, -1)， 色(0, 1, 0)
        -1, -1, 0, 0, 1 # 頂点(-1, -1)，色(0, 0, 1)
    ], dtype=GLfloat)

    # GPU上にバッファを生成
    vbo = glGenBuffers(1)
    glBindBuffer(GL_ARRAY_BUFFER, vbo)
    glBufferData(GL_ARRAY_BUFFER, data.itemsize * data.size, (GLfloat * data.size)(*data), GL_STATIC_DRAW)

    # バッファのデータとシェーダの変数を関連付け
    vao = glGenVertexArrays(1)
    glBindVertexArray(vao)
    glEnableVertexAttribArray(0)
    glEnableVertexAttribArray(1)
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 5, GLvoidp(0))
    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 5, GLvoidp(sizeof(GLfloat) * 2))
    glBindVertexArray(0)

    glClearColor(0, 0, 0, 1)
    while glfw.window_should_close(window) == glfw.FALSE:
        glClear(GL_COLOR_BUFFER_BIT)

        # 描画
        program.use()
        glBindVertexArray(vao)
        glDrawArrays(GL_TRIANGLES, 0, 3)
        glBindVertexArray(0)
        program.unuse()

        glfw.swap_buffers(window)
        glfw.wait_events()
    glDeleteVertexArrays(1, [vao])
    glDeleteBuffers(1, [vbo])
    glDeleteProgram(program.handle)

    glfw.terminate()

if __name__ == "__main__":
    main()

これを実行するためにはnumpy, PyOpenGL, glfwが必要になるので，以下のコマンドでインストールしましょう．以降のコードでも必要なライブラリはこれらのみです．

$ pip install numpy PyOpenGL glfw

　では，上のコードについて解説していきます．
　モダンなOpenGLでは以下のようなフローで描画を行います．

• GPU上にバッファを生成
  　glGenBuffersからglBufferDataまでの間では，GPU上にバッファを生成しCPU上（numpy）で構築した頂点データをGPUに送っています．これはCPUの頂点データにGPUからアクセスするためにGPUへそのデータを転送する必要があり，その転送にはかなりの時間を要するからです．描画の度に転送したり頂点毎にデータを転送するよりも，あらかじめGPUにすべてのデータを転送しておき描画の際にはGPU上にあるデータを参照する方が速いという理屈ですね．

• シェーダの生成
  　OpenGLには，プログラマブルシェーダという概念が存在します．詳細については後述しますが，簡単に言えばシェーダとは描画の仕方をプログラムによって記述したものです．プログラムによって記述することでベンダーが提供したありあわせの固定機能よりも遥かに柔軟で豊かな表現が可能になります．

• バッファのデータとシェーダの変数を関連付け
  　このフェーズでは，GPU上に配置されたデータをどのように解釈するかを決定します．三角形描画の例では一つの頂点データは5つのfloat値で構成されていて，前半2つが頂点位置，後半3つが頂点の色，これが3つ並んでいる，といった具合になります．このようなデータの構造をインターリーブと言います．ココらへんの解説は別の方の記事に任せます．

• 描画
  　このフェーズではどのバッファを使うか，どのようなシェーダ変数との関連付けでバッファのデータを使うかを指定した後，描画する頂点の数を指定して決められたプリミティブを描画します．

　このようなフェーズを経て描画されたものが以下の図になります．

プログラマブルシェーダ

　次に前述したシェーダについて説明していきます．OpenGLではバーテックスシェーダ，ジオメトリシェーダ，フラグメントシェーダ，テッセレーション制御シェーダ，テッセレーション評価シェーダ，コンピュートシェーダの6種類のシェーダがあります．この内テッセレーション制御・評価シェーダについては今回使わないため省略します．これらのシェーダの中でコンピュートシェーダ以外は描画に直接関わるシェーダであり，コンピュートシェーダは主にGPUの計算能力を活かしたデータの並列処理を行います．これらの処理を柔軟に記述するためにGLSLというC言語に似た特殊な言語を使用します．GLSLでは基本的に以下のような形に沿って記述していきます(コンピュートシェーダは少々異なります)．

#version 450 　　　　<--- 使用するGLSLのバージョンに合わせて変更する

layout (location = 0) in in_type1 in_variable1;
layout (location = 1) in in_type2 in_variable2;
             :
             :

layout (location = 0) out out_type1 out_variable1;
layout (location = 1) out out_type2 out_variable2;
             :
             :

// 関数が定義できます．
out_type1 change(in_type1 v) {
    return /* vを加工 */
}

// もちろんグローバル変数も定義できます．
const float x = 1.0;

void main() {
    out_variable1 = change(in_variable1); // 別の値に加工して渡せます．
    out_variable2 = in_variable2;         // もちろん加工する必要がなければそのまま渡せます．
    gl_* = vec4(x);                       // シェーダによってはgl_から始まる組み込みの変数があります．
}

　前述した複数のシェーダそれぞれについてGLSLで上のようなコードを書いていくわけですが，各シェーダが呼び出される順番は決まっていて，バーテックスシェーダ，ジオメトリシェーダ，フラグメントシェーダの順に呼ばれます．ここで重要なのが，outで指定された変数に値を代入することで次に呼び出されるシェーダに値を渡せるということです．例えばバーテックスシェーダでlayout (location = 0) outで指定された変数に値を代入したら，ジオメトリシェーダかフラグメントシェーダでlayout (location = 0) inで指定された変数にその値が入っているわけです(厳密に言えば，ジオメトリシェーダでは配列として受け取りますし，フラグメントシェーダでは補間が働くため違う値が入るわけですが)．
　では以上のことを踏まえて各シェーダについて説明していきましょう．

バーテックスシェーダ

　このシェーダは描画する図形の各頂点を決定します．main関数は各頂点毎に並列で呼び出されるため，それぞれの頂点をどう加工するかは同一のコードで書くことになります．またバーテックスシェーダではgl_Positionという組み込みの変数が用意されており，この変数に渡したい頂点を代入することで頂点の位置が決定されます．三角形描画で使ったバーテックスシェーダの例を示しますが，この例では頂点を加工する必要が無いためそのままの値を渡しています．

#version 450

layout (location = 0) in vec2 vertex;
layout (location = 1) in vec3 color;

layout (location = 0) out vec3 outColor;

void main() {
    gl_Position = vec4(vertex, 0, 1);
    outColor = color;
}

　この例では何も加工せずそのままの値を渡していますが，例えば三角形の位置を変えたければ平行移動行列を掛ければ良いです．ココらへんの解説は本記事の内容から逸脱するので別の解説記事に譲ります．

ジオメトリシェーダ

　このシェーダは一言で言えば，図形の変更です．例えばもともとは点を何個か描画するつもりだったけれど，もっと複雑な図形で同じものを何個か描画するということができるわけです．三角形描画の例では使う必要がなかったので出てこなかったのですが，以下のようなものです．

#version 450

layout (points) in;
layout (triangle_strip, max_vertices = 3) out;

layout (location = 0) in vec3 vertexColor[];

layout (location = 0) out vec3 outColor;

void main() {
    gl_Position = gl_in[0].gl_Position;
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = gl_in[0].gl_Position + vec4(0.1, -0.1, 0, 0);
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = gl_in[0].gl_Position + vec4(-0.1, -0.1, 0, 0);
    outColor = vertexColor[0];
    EmitVertex();
    EndPrimitive();
}

　この例では一つの点から一つの三角形に変更されるのですが，gl_Positionに値を書き込んでからEmitVertex()を呼び出すことで頂点を決定することができます．3つの頂点をすべて決めた後EndPrimitive()を呼び出すことで頂点書き込みの終了を知らせることができます．

フラグメントシェーダ

　このシェーダでは描画する図形のピクセル単位の色を決めることができます．またバーテックスシェーダやジオメトリシェーダから渡された値は特別に何らかの指定をしない限り線形補間されています．そのため，前出の三角形の図のようにグラデーションが現れるわけです．このシェーダでも前シェーダから渡された値を加工することができますが，三角形描画の例では何もせずそのまま渡しています．

#version 450

layout (location = 0) in vec3 color;

layout (location = 0) out vec4 outColor;

void main() {
    outColor = vec4(color, 1);
}

コンピュートシェーダ

　このシェーダは少し特殊です．コンピュートシェーダを理解するためにはGPUがどのように計算しているかを理解する必要があります．
　まずCPUとの違いについてですが，CPUはマルチコアプロセッサでもコアを数個から数十個程度しか持っていないのに対し，GPUは数千個持っているのが普通です．さらにCPUでは分岐予測が働き複雑な命令でも高速に実行できる一方，GPUでは複雑な処理が低速になってしまうものの数千個ものコアを活かして同一の単純な処理を並列して何万個も行うことができます．それらの処理一個一個をスレッドと呼ぶわけですが，GLSLのコンピュートシェーダではそれらをワークグループという決まった大きさの3次元のかたまりにまとめ，更にそれを複数個集めて3次元で管理する，という方法をとっています．例えば10万個の同じ処理を並列してGPUに行わせたい場合は，1個のワークグループのサイズを(32, 32, 1)とし，そのワークグループを128個集めれば良いわけです．
　以下がコンピュートシェーダの例です．この例では10万個の粒子に重力をかけて落とすということを行っています．

#version 460

// ワークグループのサイズ
layout (local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

// 粒子を表す構造体
struct particle_t {
    vec2 position;
    vec2 velocity;
};

// 10万個の粒子が格納された配列
layout (std430, binding = 0) buffer particles {
    particle_t particle[];
};

const vec2 g = vec2(0, -0.1);
const float dt = 0.1;

void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    // 実行されているスレッドのインデックス
    uint id = gl_LocalInvocationIndex + work_id * gl_WorkGroupSize.x * gl_WorkGroupSize.y * gl_WorkGroupSize.z;

    // 粒子の状態を更新
    particle[id].position += particle[id].velocity * dt;
    particle[id].velocity += g * dt;
}

　これを含めたPythonのコードを実行すると以下のGIFのようになります．[-1, 1]✕[-1, 1]の範囲に10万個の粒子があるため最初の方はカオスですが，しばらく経つと粒子が落下していく様子がはっきりと見えると思います．

　上のコードで行っていることは単純でスレッドのインデックスを求めた後，着目している粒子の位置と速度を更新しています．また粒子の状態の読み書きはシェーダーストレージバッファオブジェクトに対して行っています．シェーダーストレージバッファオブジェクトとはシェーダから自由に値の読み書きができるバッファのことです．コンピュートシェーダ，シェーダーストレージバッファオブジェクト，これらを組み合わせて使うことで非常に高速に動作するマルチエージェントシステムが作れます．ではここから本題に入りましょう．

人工生命モデル

　今回扱う人工生命モデルはGame of Life，Boids，Gray-Scottです．それぞれについて概要と結果を示していきましょう．

Game of Life

　まず有名なGame of Lifeをシミュレーションしてみましょう．このモデルについて軽く解説しておくと，Game of Lifeは2次元セルオートマトンであり，ムーア近傍8個と自身の状態，与えられたルールに基づいて状態を更新するモデルです．このルールは以下の四つで表されます．

• 誕生
  　死んだセルの周りに生存しているセルがちょうど3つあればそのセルは生きたセルに更新される．

• 生存
  　生きたセルの周りに生存しているセルが2つか2つあれば更新しない．

• 過疎
  　生きたセルの周りで生存しているセルが1つ以下ならば死んだセルに更新される．

• 過密
  　生きたセルの周りで生存しているセルが4つ以上ならば死んだセルに更新される．

　これらたった4つの単純なルールでまるで生命が生きて活動しているかのような動作がみられて非常に面白いモデルです．それではGLSLのコードを見ていきましょう．

lifegame.comp

#version 450

layout (local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

layout (std430, binding = 0) buffer cells {
    int cell[];
};

layout (std430, binding = 1) buffer temps {
    int temp[];
};

uniform ivec2 rect;

int get(int x, int y) {
    return cell[x + y * rect.x];
}

int neighbor(int x, int y, int dx, int dy) {
    int _x = x + dx < 0 ? rect.x - 1 : x + dx >= rect.x ? 0 : x + dx;
    int _y = y + dy < 0 ? rect.y - 1 : y + dy >= rect.y ? 0 : y + dy;
    return get(_x, _y);
}

void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    uint id = gl_LocalInvocationIndex + work_id * 32 * 32;

    int x = int(mod(id, rect.x));
    int y = int(id / rect.x);

    int num = 0;
    num += neighbor(x, y, -1, 0);
    num += neighbor(x, y, -1, -1);
    num += neighbor(x, y, 0, -1);
    num += neighbor(x, y, 1, -1);
    num += neighbor(x, y, 1, 0);
    num += neighbor(x, y, 1, 1);
    num += neighbor(x, y, 0, 1);
    num += neighbor(x, y, -1, 1);
    temp[id] = num == 3 ? 1 : num == 2 ? cell[id] : 0;
}

copy.comp

#version 450

layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

layout (std430, binding = 0) buffer cells {
    int cell[];
};

layout (std430, binding = 1) buffer temps {
    int temp[];
};

void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    uint id = gl_LocalInvocationIndex + work_id * 32 * 32;

    cell[id] = temp[id];
}

　Game of Lifeでは二つのコンピュートシェーダを使用します．1つ目はlifegame.compでセルの更新を行っていますが，出力先は実際に表示されるcellではなくtempというバッファです．2つ目のcopy.compではtempの値をcellに代入しているだけです．呼び出す順番はlifegame.comp→copy.compです．なぜこのような冗長なことをしなければならないのでしょうか？
　これに関わるキーワードとしてデータ競合というものがあります．Game of Lifeではセルの更新に近傍の値を使います．しかしGPUでは注目しているセルの近傍を読み出すのと同時に別のスレッドでその近傍セルが更新される，ということが起こります．これがまさにデータ競合です．これが起こった時セルの値を正しく読み取ることができません．これを解決する方法としてスレッドの同期があります．しかしGLSLではワークグループ内のスレッド同期命令はあるものの，残念ながらワークグループ間でのスレッド同期命令はないようです．そのためコンピュートシェーダを二つ用意し，1つは書き込み，1つは読み込みと分けることで同じメモリに対する複数のアクセスを排除しました．
　このようにしてある時刻での各セルの状態は計算できたわけですが，次にセルを表示しなければなりません．ここからはcellに書き込まれたデータをバーテックスシェーダで取り出して各セルを正方形として描画していく，ということを行っていくのですが，困ったことにバーテックスシェーダではこの「正方形として描画」する，ということができません（厳密にはできるのですがオフセットがどうたらと非常にめんどくさいです）．そこで一旦バーテックスシェーダでは，格子点上に「点を打つ」ということだけを行いましょう．

lifegame.vert

#version 450

layout (std430, binding = 0) buffer cells {
    int cell[];
};

layout (location = 0) out vec4 vertexColor;

uniform ivec2 rect;

void main() {
    float fid = gl_VertexID;
    float x = mod(fid, rect.x) / rect.x;
    float y = int(fid / rect.x) / float(rect.y);

    gl_Position = vec4(vec2(x, y) * 2 - 1, 0, 1);
    vertexColor = vec4(cell[gl_VertexID]);
}

　ここで登場するのがジオメトリシェーダです．ジオメトリシェーダの役割を思い出すと「描画する図形の変更」です．これはまさに今必要としている「点から正方形への変換」という機能と合致します．ではジオメトリシェーダを見ていきましょう．

lifegame.geom

#version 450

layout (points) in;
layout (triangle_strip, max_vertices = 4) out;

layout (location = 0) in vec4 vertexColor[];
layout (location = 0) out vec4 outColor;

uniform ivec2 rect;

void main(void) {
    gl_Position = gl_in[0].gl_Position;
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = gl_in[0].gl_Position + vec4(2.0 / rect.x, 0, 0, 0);
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = gl_in[0].gl_Position + vec4(0, 2.0 / rect.y, 0, 0);
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = gl_in[0].gl_Position + vec4(2.0 / rect.x, 2.0 / rect.y, 0, 0);
    outColor = vertexColor[0];
    EmitVertex();
    EndPrimitive();
}

　ジオメトリシェーダで複雑な図形に変換する場合，たいていtriangle_stripに変換します．今回もtriangle_stripに変換していて，描画したい正方形の頂点を4つ入力しているだけです．以上のようにして1000x1000の2次元セルオートマトンを実行した結果が以下の図です．

Boids

　続いてBoidsモデルです．Boidsは「鳥」を意味するBirdと「もどき」を意味する接尾辞である-oidを合わせて鳥もどきという意味です．鳥や魚の群れをシミュレーションするマルチエージェントシステムであり，発表されたのが1980年代とかなり昔ですが現在でも映画などで使われるほど実際の群れの挙動をよく表したモデルです．
　このモデルでは以下の3つのルールに基づいて各エージェントが動きます．

• 集合
  　群れの中心方向に向かう力が働きます．

• 分離
  　近すぎる個体から離れる方向に力が働きます．

• 整列
  　群れ全体の向いている方向に合わせる力が働きます．

　これら3つのルールに関して，影響を及ぼす距離と力の働く割合を変化させることで群れのマクロな動きを変化させることができます．ではコンピュートシェーダのコードを示していきましょう．

boids.comp

#version 450

layout (local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

struct agent_t {
    vec2 position;
    vec2 velocity;
};

layout (std430, binding = 0) buffer agents {
    agent_t agent[];
};

layout (std430, binding = 1) buffer temps {
    agent_t temp[];
};

uniform int num_agents;
uniform float min_velocity;
uniform float max_velocity;
uniform vec2 mouse_position;

const float pi = 3.1415926535;

const float fc = 0.005;
const float fs = 0.5;
const float fa = 0.01;
const float fb = 0.04;
const float fw = 2.8;
const float dc = 0.8;
const float ds = 0.03;
const float da = 0.5;
const float dw = 0.5;
const float ac = pi / 2;
const float as = pi / 2;
const float aa = pi / 3;


void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    uint id = gl_LocalInvocationIndex + work_id * 32 * 32;

    int num_c = 0;
    int num_s = 0;
    int num_a = 0;
    vec2 f_c = vec2(0);
    vec2 f_s = vec2(0);
    vec2 f_a = vec2(0);
    for(int i = 0; i < num_agents; i++) {
        float l = length(agent[i].position - agent[id].position);
        float a = acos(dot(normalize(agent[i].position), normalize(agent[id].position)));
        f_c += l < dc && a < ac ? (num_c++, agent[i].position) : vec2(0);
        f_s += l < ds && a < as ? (num_s++, agent[id].position - agent[i].position) : vec2(0);
        f_a += l < da && a < aa ? (num_a++, agent[i].velocity) : vec2(0);
    }
    f_c = num_c == 0 ? vec2(0) : (f_c / num_c - agent[id].position) * fc;
    f_s *= fs;
    f_a = num_a == 0 ? vec2(0) : (f_a / num_a - agent[id].velocity) * fa;
    float dist_center = length(agent[id].position);
    vec2 f_b = dist_center > 1 ? -agent[id].position * (dist_center - 1) / dist_center * fb : vec2(0);

    vec2 f_w = length(agent[id].position - mouse_position) < dw ? (agent[id].position - mouse_position) * fw : vec2(0);

    temp[id].velocity = agent[id].velocity + f_c + f_s + f_a + f_b + f_w;
    float lv = length(temp[id].velocity);
    vec2 vdirection = normalize(temp[id].velocity);
    temp[id].velocity = lv < min_velocity ? vdirection * min_velocity : (lv > max_velocity ? vdirection * max_velocity : temp[id].velocity);
    temp[id].position = agent[id].position + temp[id].velocity;
}

copy.comp

#version 450

layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

struct agent_t {
    vec2 position;
    vec2 velocity;
};

layout (std430, binding = 0) buffer agents {
    agent_t agent[];
};

layout (std430, binding = 1) buffer temps {
    agent_t temp[];
};

void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    uint id = gl_LocalInvocationIndex + work_id * 32 * 32;

    agent[id] = temp[id];
}

　今回のコードでも2つのコンピュートシェーダを使っていて，役割はLife of Gameのものと全く同じです．boids.compの方についてですが，やっていることは単純で注目しているエージェントの近傍個体を調べて，それぞれの間で働く力を計算し，働く力をすべて足し合わせてから位置と速度の更新を行っています．
　次にジオメトリシェーダです．今回も点から三角形へ変換するので使用します．

boids.geom

#version 450

layout (points) in;
layout (triangle_strip, max_vertices = 3) out;

layout (location = 0) in vec4 vertexColor[];
layout (location = 1) in vec2 direction[];
layout (location = 0) out vec4 outColor;

float cross2d(vec2 a, vec2 b) 
{
    return a.x * b.y - a.y * b.x;
}

mat2 rot2d(float cs, float sn) {
    return mat2(
        cs, -sn,
        sn, cs
    );
}

void main(void) {
    float cs = dot(direction[0], vec2(1, 0));
    float sn = cross2d(direction[0], vec2(1, 0));
    vec2 position = gl_in[0].gl_Position.xy;
    mat2 rotm = rot2d(cs, sn);

    gl_Position = vec4(rotm * vec2(0, 0.004) + position, 0, 1);
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = vec4(rotm * vec2(0, -0.004) + position, 0, 1);
    outColor = vertexColor[0];
    EmitVertex();
    gl_Position = vec4(rotm * vec2(0.02, 0) + position, 0, 1);
    outColor = vertexColor[0];
    EmitVertex();
    EndPrimitive();
}

　このコードでは内積と外積を使いcosとsinを求め，その値に基づき回転行列を生成して，エージェントを表す三角形を回転させています．ジオメトリシェーダであればこのような加工をすることもできます．以上のコードで1000個体のシミュレーションを行った結果が以下の図です．マウスの入力も受け付けていて，ドラッグした位置から近くにいたエージェントが離れていきます．

Gray-Scott

　最後にGray-Scottです．反応拡散型の方程式系でセルの持つ値が二値ではなく実数値な点，ノイマン近傍である点，セルの更新ルール以外はGame of Lifeと同じです．コンピュートシェーダは以下のようになります．

gray-scott.py

#version 450

layout (local_size_x = 32, local_size_y = 32, local_size_z = 1) in;

layout (std430, binding = 0) buffer us {
    float u[];
};

layout (std430, binding = 1) buffer u_temps {
    float u_temp[];
};

layout (std430, binding = 2) buffer vs {
    float v[];
};

layout (std430, binding = 3) buffer v_temps {
    float v_temp[];
};

uniform ivec2 rect;

float get_u(int x, int y) {
    return u[x + y * rect.x];
}

float neighbor_u(int x, int y, int dx, int dy) {
    int _x = x + dx < 0 ? rect.x - 1 : x + dx >= rect.x ? 0 : x + dx;
    int _y = y + dy < 0 ? rect.y - 1 : y + dy >= rect.y ? 0 : y + dy;
    return get_u(_x, _y);
}

float get_v(int x, int y) {
    return v[x + y * rect.x];
}

float neighbor_v(int x, int y, int dx, int dy) {
    int _x = x + dx < 0 ? rect.x - 1 : x + dx >= rect.x ? 0 : x + dx;
    int _y = y + dy < 0 ? rect.y - 1 : y + dy >= rect.y ? 0 : y + dy;
    return get_v(_x, _y);
}

const float dx = 0.01;
const float dt = 1;
const float Du = 2e-5;
const float Dv = 1e-5;
const float f = 0.022;
const float k = 0.051;

void main() {
    uint work_id = gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y + 
                   gl_WorkGroupID.y * gl_NumWorkGroups.x +
                   gl_WorkGroupID.x;
    uint id = gl_LocalInvocationIndex + work_id * 32 * 32;

    int x = int(mod(id, rect.x));
    int y = int(id / rect.x);

    float laplacian_u = (neighbor_u(x, y, 1, 0) + neighbor_u(x, y, -1, 0) + neighbor_u(x, y, 0, 1) + neighbor_u(x, y, 0, -1) - 4 * u[id]) / (dx * dx);
    float laplacian_v = (neighbor_v(x, y, 1, 0) + neighbor_v(x, y, -1, 0) + neighbor_v(x, y, 0, 1) + neighbor_v(x, y, 0, -1) - 4 * v[id]) / (dx * dx);
    float dudt = Du * laplacian_u - u[id] * v[id] * v[id] + f * (1 - u[id]);
    float dvdt = Dv * laplacian_v + u[id] * v[id] * v[id] - (f + k) * v[id];
    u_temp[id] = u[id] + dt * dudt;
    v_temp[id] = v[id] + dt * dvdt;
}

　copy.compはGame of Lifeのものとほぼ同じなので省略します．Gray-Scottは2チャネルで構成された2次元セルオートマトンでチャネル間の相互作用により複雑な挙動が生まれるモデルであるため，登場するバッファもu, vの2種類あります．
　バーテックスシェーダ，ジオメトリシェーダについてもGame of Lifeと同一のものなので省略します．以上のようにして描画した結果は以下のようになります．

　このモデルの挙動は非常に面白いので，最終的にどうなるかはご自身で確かめてみてください．

おわりに

　本記事では，OpenGLを使ったALifeの高速なシミュレーションを実装し紹介しました．今回お見せしたGIF画像はすべてCPU内臓のGPUを使ってシミュレーションしたものであり，強力なGPUがなくともこれだけ速く動作します．興味のある方は，今回登場しなかった人工生命モデルをご自身で実装してみるのも良いかもしれません．比較的新しい環境であれば今回紹介した技術は大体使えると思います．