【渋滞学】セル・オートマトンの実装

頭痛のひどいRockinWoolです。しかし、時期的にはかなり忙しいのでプログラミングもはかどっています。最近の趣味は公園で散歩をしながら肩をほぐすことですね。
さて、ちょっと前に渋滞学という本を読んだのですが、その中にあったセル・オートマトンの実装をずっとやってみたいなと思いながらやっていなかったのでやりました。今回はpybind11を使ってセル・オートマトン部分をC++で作り、作図部分やデータ整理部分をpythonに任せる形にしています。実装はGithubにこっそり上げて置きます。

プログラムの構成

まず、メインプログラムは実験を行ってデータをcsvファイルに出力するmain.pyとcsvを元に作図をするdraw.pyの2つから構成することにしました。このようにすれば、main.pyはpybind11を使って連携するがdraw.pyはピュアpythonの実装で良くなるので便利でした。ざっとフォルダ階層を書くと下のような感じです。

-main.py
-draw.py
-CMakeLists.txt
-CA.cpp
-CA.hpp
-utils
  |-caHistory.cpp

CA.cppはCA.hppで定義された関数をpythonから使用できるようにしていて、細かい実装はcaHistory.cppに書いてあります。余談ですが、CAはCellAutomatonの略、historyとつけているのはセル・オートマトンでどの位置にエージェントがいたかどうかを返すことを明確にしたかったからです。

main.pyの概要

まずは基本となるmain.pyの解説から。難しい処理はC++に任せてしまっているので中身は超単純です。1~100個のエージェントを召喚してセル・オートマトンをして、その情報を保存することを2回繰り返しています。1回目は1個から始めて100個まで増やして混雑具合を見ています。2回目はその逆で100個から減らしていってます。cl.caHistory()でセル・オートマトンの実験をしていますが、その時に計算される混雑の回数congestionはゲッターであるget_congestion()を実装して入手する形になっています。これはpybind11の仕様でメンバ変数を直接参照できないことと、メンバ変数は基本的にprivateにしてゲッターを使ってアクセスするというベストプラクティスに沿った実装をしていることの2点があります。

Number_of_cells =100
Number_of_trial =1000
congestions = []
cl = ca.CA(Number_of_cells)
for i in range(Number_of_cells):
    ret = cl.caHistory(Number_of_trial)
    congestions.append(cl.get_congestion())
    cl.add_agent()

for i in range(Number_of_cells):
    cl.remove_agent()
    ret = cl.caHistory(Number_of_trial)
    congestions.append(cl.get_congestion())

df = pd.DataFrame(congestions)
df.to_csv("./data/congestions.csv")

それでは、このmain.pyで使用されるcaHistory(), get_congestion(), add_agent(), remove_agent()を実装していきましょう。

CA.cppの実装

こちらも中身は単純で、pybindのincludeを行う部分と関数をpython側に共有させる部分の2つから構成されています。書き方が少し独特ですが4つの関数が公開されるのがわかります。

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
#include <pybind11/complex.h>
#include <pybind11/functional.h>
#include <pybind11/chrono.h>
#include "CA.hpp"
namespace py = pybind11;

PYBIND11_MODULE(cellauto, m)
{
    py::class_<CA>(m, "CA")
        .def(py::init<int>())
        .def("caHistory", &CA::caHistory)
        .def("add_agent", &CA::add_agent)
        .def("remove_agent", &CA::remove_agent)
        .def("get_congestion", &CA::get_congestion);
}

CA.hppの内容

先程の4つの関数のプロトタイプ宣言を行っている部分です。ちなみに前回の記事で書きましたとおり、コンストラクタも作らないといけないのでpublicには5つの関数が登録されています。またprivateにはそれらの関数が使う情報を登録しています。update()関数はセル・オートマトンの具体的動作を規定していて、前のセルにエージェントがいたら止まるなどの挙動を定義しています。これをcaHistory()ではint回呼び出してシミュレートするわけですね。

#include <iostream>
#include <string>
#include <vector>
using namespace std;

class CA{
    public:
    CA(int);
    vector<vector<int>> caHistory(int);
    void add_agent();
    void remove_agent();
    vector<int> get_congestion();

    private:
    int size;
    vector<int> cells;
    vector<int> congestion;
    vector<vector<int>> history;
    void update();
};

caHistory.cppの実装

記事のほとんどがこのプログラムの内容で埋まってしまっていますが、実際はかなり単純なプログラムです。ほとんどの分量はupdate()関数の条件文で埋まっているのでcaHistory(), add_agent(), remove_agent(), get_congestion()の本文は2~5行程度になっています。それぞれ与えられた役割が動作するようになっています。

#include "../CA.hpp"

CA::CA(int size):size(size){
    this->cells.resize(this->size, 0);
    this->congestion.resize(this->size,0);
    this->cells[0]=1;
    this->history.push_back(cells);
}

void CA::update() {
    std::vector<int> newCells(this->size, 0);
    for (int i = 1; i < this->size - 1; ++i) {
        // 簡単なルール：前後のセルの状態に基づいて次の状態を更新
        if      (this->cells[i+1]==1 && this->cells[i]==1 && this->cells[i-1]==1){newCells[i]=1;this->congestion[i]++;}
        else if (this->cells[i+1]==1 && this->cells[i]==1 && this->cells[i-1]==0){newCells[i]=1;this->congestion[i]++;}
        else if (this->cells[i+1]==1 && this->cells[i]==0 && this->cells[i-1]==1){newCells[i]=1;}
        else if (this->cells[i+1]==1 && this->cells[i]==0 && this->cells[i-1]==0){newCells[i]=0;}
        else if (this->cells[i+1]==0 && this->cells[i]==1 && this->cells[i-1]==1){newCells[i]=0;}
        else if (this->cells[i+1]==0 && this->cells[i]==1 && this->cells[i-1]==0){newCells[i]=0;}
        else if (this->cells[i+1]==0 && this->cells[i]==0 && this->cells[i-1]==1){newCells[i]=1;}
        else if (this->cells[i+1]==0 && this->cells[i]==0 && this->cells[i-1]==0){newCells[i]=0;}
        else{cout<<"error"<<endl;}
    }
    // TODO; 両側の境界条件を書きましょう----//
    if      (this->cells[1]==1 && this->cells[0]==1 && this->cells[this->size-1]==1){newCells[0]=1;this->congestion[0]++;}
    else if (this->cells[1]==1 && this->cells[0]==1 && this->cells[this->size-1]==0){newCells[0]=1;this->congestion[0]++;}
    else if (this->cells[1]==1 && this->cells[0]==0 && this->cells[this->size-1]==1){newCells[0]=1;}
    else if (this->cells[1]==1 && this->cells[0]==0 && this->cells[this->size-1]==0){newCells[0]=0;}
    else if (this->cells[1]==0 && this->cells[0]==1 && this->cells[this->size-1]==1){newCells[0]=0;}
    else if (this->cells[1]==0 && this->cells[0]==1 && this->cells[this->size-1]==0){newCells[0]=0;}
    else if (this->cells[1]==0 && this->cells[0]==0 && this->cells[this->size-1]==1){newCells[0]=1;}
    else if (this->cells[1]==0 && this->cells[0]==0 && this->cells[this->size-1]==0){newCells[0]=0;}
    else{cout<<"error"<<endl;}
    if      (this->cells[0]==1 && this->cells[this->size-1]==1 && this->cells[this->size-2]==1){newCells[this->size-1]=1;this->congestion[this->size-1]++;}
    else if (this->cells[0]==1 && this->cells[this->size-1]==1 && this->cells[this->size-2]==0){newCells[this->size-1]=1;this->congestion[this->size-1]++;}
    else if (this->cells[0]==1 && this->cells[this->size-1]==0 && this->cells[this->size-2]==1){newCells[this->size-1]=1;}
    else if (this->cells[0]==1 && this->cells[this->size-1]==0 && this->cells[this->size-2]==0){newCells[this->size-1]=0;}
    else if (this->cells[0]==0 && this->cells[this->size-1]==1 && this->cells[this->size-2]==1){newCells[this->size-1]=0;}
    else if (this->cells[0]==0 && this->cells[this->size-1]==1 && this->cells[this->size-2]==0){newCells[this->size-1]=0;}
    else if (this->cells[0]==0 && this->cells[this->size-1]==0 && this->cells[this->size-2]==1){newCells[this->size-1]=1;}
    else if (this->cells[0]==0 && this->cells[this->size-1]==0 && this->cells[this->size-2]==0){newCells[this->size-1]=0;}
    else{cout<<"error"<<endl;}
    // --------------------------------------//
    this->cells = newCells;
    this->history.push_back(this->cells);
}

vector<vector<int>> CA::caHistory(int generation){
    // 10世代分進化させながら表示
    for (int gen = 0; gen < generation; ++gen) {
        update();
    }
    return this->history;
}

void CA::add_agent(){
    // clear the information of congestion before add agent.
    this->congestion = vector<int>(this->size,0);
    for(int i=0;i<this->size;i++){
        if(this->cells[i]==0){
            this->cells[i]=1;
            break;
        }
    }
}

void CA::remove_agent(){
    // clear the information of congestion before add agent.
    this->congestion = vector<int>(this->size,0);
    for(int i=0;i<this->size;i++){
        if(this->cells[i]==1){
            this->cells[i]=0;
            break;
        }
    }
}

vector<int> CA::get_congestion(){
    return this->congestion;
}

draw.pyによる描画結果

draw.pyの実装は大した内容では無いので割愛します。main.pyを実行して得られたcsvファイルをベースにdraw.pyが書いた図が下記です。

まとめ

今回はpybind11を使ってセル・オートマトンを実装してみました。内容が少しヘビーになってしまったかもしれませんが、興味のある方はコードの改善点などドシドシ送ってください。それでは。