JavaScriptでGray-Scottモデルを実装してみた話

はじめに

最近、人工生命に関する本を読んだのですが、その中で「Gray-Scottモデル」という面白いものがあったのでJavaScriptとCanvasで作ってみました。WebGLを使っていないので描画処理がちょっと重いかも。。。

コチラです。→Gray-Scottモデル v1

追記(2020/04/04):WebGL2を使って改良しました。↓

WebGL2でGray-Scottモデルを改良してみた

うまくいくとこんな縞模様や...

縞模様のチューリングパターン

増殖する斑点模様など...

斑点模様のチューリングパターン

様々な模様が自然に形成されていくのを見ることができます。これらの模様はチューリングパターンというそうで、シマウマや熱帯魚、貝などの自然界に存在する生物においても確認することができます。

参考にした本はコチラです。↓

作って動かすALife ―実装を通した人工生命モデル理論入門 (Amazon)

以下の説明やコードはこの本を参考にしています。Gray-Scottモデル以外にも、ニューラルネットワークや遺伝的アルゴリズムを使った人工生命へのアプローチが載っているので面白いです。

反応拡散方程式

その昔、アラン・チューリングという天才はこのような模様の形成を物質の反応と拡散を表す方程式によって数学的に表せることを示しました。これを「反応拡散系」と呼び、今回作ったGray-Scottモデルはその一種です。物質Uと物質Vの濃度\(u, v\)の偏微分方程式として表されます。

\[ \frac{\partial u}{\partial t} = D_u \nabla^2 u - uv^2 + f(1-u) \] \[ \frac{\partial v}{\partial t} = D_v \nabla^2 v + uv^2 - v(f+k) \]

第1項が物質の拡散(\(D_u > D_v\)のため物質Uのほうが速く拡散する)、第2項が「U + 2V → 3V」という物質の反応を表し、最後の項がUの流入(あるいはVの流出)を表しています。また、定数\(f, k\)の値を変化させることで生成されるパターンの種類を変化させることができます。先述した自作Gray-Scottモデルでも、それぞれFeed, Killというパラメータとして変化できるようにしました。

JavaScriptで実装

参考書籍ではPythonとVisPyというライブラリで実装されていたのですが、今回はJavaScriptとCanvasを使って実装してみました(この辺の事情は後述)。

gray_scott.js
let interval_id1 = 0;
let interval_id2 = 0;
const GRID_X = 150;         // 参考書籍では256x256
const GRID_Y = 150;
const SQUARE_SIZE = 10;
const dx = 0.01707;         // 参考書籍では0.01なので(256/150)倍してみた
const dt = 3;               // 参考書籍では1だがスピードアップのため3に
const Du = 2e-5;
const Dv = 1e-5;
const presets = [
    { feed: 0.022, kill: 0.051 },   // stripe                               
    { feed: 0.035, kill: 0.065 },   // spots
    { feed: 0.012, kill: 0.050 },   // wandering bubbles
    { feed: 0.025, kill: 0.050 },   // waves
    { feed: 0.040, kill: 0.060 },   // amorphous
    // 以上5つは参考書籍より。以下は勝手に追加。
    { feed: 0.025, kill: 0.060 },   // waving spots
    { feed: 0.030, kill: 0.060 },   // snapping strings
    { feed: 0.011, kill: 0.046 }    // balloons
]
let f = presets[0].feed;
let k = presets[0].kill;
let u = Array(GRID_Y);
let v = Array(GRID_Y);
let lap_u = Array(GRID_Y);
let lap_v = Array(GRID_Y);
for (let i = 0; i < GRID_Y; i++) {
    u[i] = Array(GRID_X);
    v[i] = Array(GRID_X);
    lap_u[i] = Array(GRID_X);
    lap_v[i] = Array(GRID_X);
}

// JQuery UIでUIを実装してます
$(function() {
    $("#init_btn").button().click(function () {
        if (interval_id1 == 0) {
            $("#ctrl_btn").text("停止");
        } else {
            clearInterval(interval_id1);
            clearInterval(interval_id2);
            interval_id1 = 0;
            interval_id2 = 0;
        }
        init();
    });
    $("#init_btn").text("初期化");
    $("#ctrl_btn").button().click(function (event) {
        if (interval_id1 == 0) {
            interval_id1 = setInterval(update, 4);
            interval_id2 = setInterval(draw, 100, u);
            $("#ctrl_btn").text("停止");
        } else {
            clearInterval(interval_id1);
            clearInterval(interval_id2);
            interval_id1 = 0;
            interval_id2 = 0;
            $("#ctrl_btn").text("再開");
        }
    });
    $("#ctrl_btn").text("停止");
    $("#preset_select").change(function (event) {
        $("#feed_slider").slider("value", presets[event.target.value].feed);
        $("#kill_slider").slider("value", presets[event.target.value].kill);
    });
    $("#feed_slider").slider({
        value: f, min: 0.000, max: 0.050, step: 0.001,
        change: function (event, ui) {f = ui.value; $("#feed_span").text(f.toFixed(3));},
        slide: function (event, ui) {f = ui.value; $("#feed_span").text(f.toFixed(3));}
    });
    $("#feed_slider").slider("value", f);
    $("#kill_slider").slider({
        value: k, min: 0.040, max: 0.075, step: 0.001,
        change: function (event, ui) {k = ui.value; $("#kill_span").text(k.toFixed(3));},
        slide: function (event, ui) {k = ui.value; $("#kill_span").text(k.toFixed(3));}
    });
    $("#kill_slider").slider("value", k);

    init();
});

function init() {
    // initialize u, v
    for (let i = 0; i < GRID_Y; i++) {
        u[i].fill(1);
        v[i].fill(0);
    }

    const x_0 = Math.floor(GRID_X / 2) - Math.floor(SQUARE_SIZE / 2);
    const y_0 = Math.floor(GRID_Y / 2) - Math.floor(SQUARE_SIZE / 2);
    for (let i = 0; i < SQUARE_SIZE; i++) {
        u[y_0 + i].fill(0.5, x_0, x_0 + SQUARE_SIZE);
        v[y_0 + i].fill(0.25, x_0, x_0 + SQUARE_SIZE);
    }

    for (let i = 0; i < GRID_Y; i++) {
        for (let j = 0; j < GRID_X; j++) {
            u[i][j] += Math.random()*0.09;
            v[i][j] += Math.random()*0.09;
        }
    }

    draw(u);
    interval_id1 = setInterval(update, 4);
    // draw処理は負荷が高いので間隔を長めにしました
    interval_id2 = setInterval(draw, 125, u);
}

function update() {
    // calculate laplacian
    for (let i = 0; i < GRID_Y; i++) {
        for (let j = 0; j < GRID_X; j++) {
            const dec_x = (j - 1 < 0) ? GRID_X - 1 : j - 1;
            const inc_x = (j + 1 > GRID_X - 1) ? 0 : j + 1;
            const dec_y = (i - 1 < 0) ? GRID_Y - 1 : i - 1;
            const inc_y = (i + 1 > GRID_Y - 1) ? 0 : i + 1;
            lap_u[i][j] = (u[dec_y][j] + u[inc_y][j] + u[i][dec_x] + u[i][inc_x] - 4*u[i][j]) / (dx*dx);
            lap_v[i][j] = (v[dec_y][j] + v[inc_y][j] + v[i][dec_x] + v[i][inc_x] - 4*v[i][j]) / (dx*dx);
        }
    }

    // update u, v
    for (let i = 0; i < GRID_Y; i++) {
        for (let j = 0; j < GRID_X; j++) {
            let du = Du*lap_u[i][j] - u[i][j]*v[i][j]*v[i][j] + f*(1 - u[i][j]);
            let dv = Dv*lap_v[i][j] + u[i][j]*v[i][j]*v[i][j] - v[i][j]*(f + k);
            u[i][j] += du*dt;
            v[i][j] += dv*dt;
        }
    }
}

function draw(mat) {
    const canvas = document.getElementById('canvas');
    const cell_width = Math.floor(canvas.width / GRID_X);
    const cell_height = Math.floor(canvas.height / GRID_Y);

    if (canvas.getContext) {
        const ctx = canvas.getContext('2d');
        for (let i = 0; i < GRID_Y; i++) {
            for (let j = 0; j < GRID_X; j++) {
                const x = Math.floor(cell_width * j);
                const y = Math.floor(cell_height * i);

                let val = mat[i][j];
                if (val < 0) {
                    val = 0;
                    mat[i][j] = 0;
                } else if (val > 1) {
                    val = 1;
                    mat[i][j] = 1;
                }
                const l = Math.floor(val * 100);

                ctx.fillStyle = 'hsl(0, 0%,' + l + '%)';
                ctx.fillRect(x, y, cell_width, cell_height);
            }            
        }
}
}

おわりに

今回PythonではなくJavaScriptを用いたのは、ネットに公開するのが容易であることももちろんですが一番最初の理由はOpenGLの知識を使わずに済むと思ったからです。参考書籍で使われているVisPyはOpenGLの知識が必要だし、他の可視化ライブラリであるmatplotlibやbokehはリアルタイムでプロットするには向いていない。。。というわけでJavaScriptとCanvasを使って可視化することにしました。しかし、今回実装してみてわかったのが、このままでは描画処理の負荷がかなり高いためWebGLを使った2D描画が必要だということです。結局OpenGLは避けられないのか。。。

ちなみにWebGLを使ってGray-Scottモデルを可視化しているサイトもありました。↓

WebGL Gray-Scott Explorer

というわけで、今後WebGLを使った修正版を作るかもしれません。また、他にも人工生命関連のことをいろいろ試していきたいです。

追記(2020/04/04):改良しました。↓

WebGL2でGray-Scottモデルを改良してみた