SimPyで離散事象シミュレーション(4) アニメーションを追加しよう

# はじめに

SimPyというPythonの離散事象シミュレーション用のパッケージを見つけて試してみたら気に入ったので自分用の備忘録も兼ねて使い方をまとめていく．

* [SimPyのドキュメント](https://simpy.readthedocs.io/en/latest/)
* [SimPyのソース](https://bitbucket.org/simpy/simpy/src/default/)
* [第1回の記事](https://qiita.com/j54854/items/687f8842be1d3e212384)（はじめの一歩）
* [第2回の記事](https://qiita.com/j54854/items/816d7b974562fc4e93bb)（リソースを理解しよう：Container編）
* [第3回の記事](https://qiita.com/j54854/items/2bcc189acb0cd9667857)（リソースを理解しよう：StoreとResource編）

第4回目の今回は，シミュレーションの途中経過のアニメーションに挑戦しよう．といっても，SimPy自体にはアニメーションの機能は存在しないので別の手段の助けを借りることになる．Pythonにはさまざまなビジュアライゼーションのツールが存在するので，基本的には好みのものを利用すればよい．ここでは，それらのうち，

* [matplotlib](https://matplotlib.org/)
* [p5（p5py）](https://p5.readthedocs.io/en/latest/)

の2つの手段を用いた例をそれぞれ簡単に紹介することにする．

# matplotlibを用いた例

第2回に取り上げた在庫管理モデルの例を思い出そう．Containerリソースをデフォルトのまま用いた例である．リソースの属性`level`を参照することで在庫量がわかるので，その推移を折れ線グラフで表示してみる．

```python
import numpy as np
import matplotlib.pyplot as plt
import random
import simpy

def manager(env):
    env.model.ordered = False  # no back order to receive
    env.stocktake = env.event()  # create the first signal (event)
    while True:
        yield env.stocktake
        report(env)
        if not env.model.ordered and env.model.level <= 10:
        # only when no back order to receive
        # reorder point = 10
            env.process(deliverer(env))  # activate deliverer
            env.model.ordered = True  # back order will be received
        env.stocktake = env.event()  # create the next signal (event)

def deliverer(env):
    yield env.timeout(5)  # delivery lead time = 5
    env.model.put(20)  # back order is recieved
    env.model.ordered = False  # no back order to receive

def customer(env):
    while True:
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        how_many = random.randint(1, 3)
        env.model.get(how_many)
        env.stocktake.succeed()  # signal for stocktaking (event)

def report(env):
    print('[{}] current level: {}, orderd: {}, queue length: {} '.format(round(env.now), env.model.level, env.model.ordered, len(env.model.get_queue)))

# ---------- code for visualization ----------
    env.x.append(env.now)  # time
    env.y.append(env.model.level)  # container's level
# ---------- ---------- ---------- ----------

def main():
    env = simpy.Environment()
    env.model = simpy.Container(env, init=10)  # model is marely a Container
    env.process(manager(env))
    env.process(customer(env))
#    env.run(until=200)

# ---------- code for visualization ----------
    env.x = [0]  # time
    env.y = [10]  # container's level
    fig = plt.figure(1)
    ax = fig.add_subplot(111)
    ax.set_xlabel('time')
    ax.set_ylabel('number')
    ax.set_xlim(0, 200)
    ax.set_ylim(0, 30)
    line, = ax.plot(env.x, env.y, label='level')
    ax.legend()
    ax.grid()
    for t in range(1, 200):
        env.run(until=t)  # stepwise execution
        line.set_data(env.x, env.y)
        plt.pause(0.1)
    plt.show()
# ---------- ---------- ---------- ----------

if __name__ == "__main__":
    main()
```

`main()`関数の後半に折れ線グラフを表示するコードが追加されている．ポイントは，for文で時刻tを1から200まで1ずつ増加させ，各ループで`env.run(until=t)`を呼んでいることである．これによって時間1刻みでシミュレーションを進め，そのたびに折れ線グラフを更新していっている．

なお，表示に使うデータ（`env.x`，`env.y`）は，このfor文の中で更新してもいいが，時間1よりも細かい単位での変化を反映させるために，上の例では，関数`report()`の中で更新している．

# p5（p5py）を用いた例

上の例のように，折れ線グラフなどのグラフ表示でシミュレーションの経緯を追いかけたい場合は，matplotlibが便利である．しかし，ときには，もう少し自由度の高いアニメーションを作成したくなることもある．その場合は，いわゆるドロー系のツールに頼りたい．

ここでは，Pythonで利用できるドロー系のパッケージの一例として，p5（p5py）を利用してみる．これは，[Processing](https://processing.org/)をPythonに移植したものだと考えればいいだろう．もしかするとあまりメジャーなパッケージではないかもしれないが，（筆者のような）Processingユーザにとってはありがたい選択肢である．

p5（p5py）とProcessingの違いは[ここ](https://p5.readthedocs.io/en/latest/guides/for-processing-users.html)にまとめられている．なお，余談として，Processing自体にも[Python Mode](https://py.processing.org/)とよばれる機能があるが，これはJythonなので，SimPyと組み合わせるのは難しい．

p5（p5py）の基本的な構文は下のようになる．

```python
from p5 import *

def setup():
    size(**スケッチの幅**, **スケッチの高さ**)

# setup()関数の中に初期設定コードを書く


def draw():
    background(**背景色の指定**)

# draw()関数の中に描画のコードを書く


if __name__ == '__main__':
    run()  # 実際に描画する際にrun（）関数を呼ぶ
```

`draw()`関数はループ実行されるので，これをアニメーションのフレームに対応させる．最後に，描画のためのp5pyの`run()`関数を呼ぶ必要があるので，SimPyの方の`run()`関数は`draw()`関数の中でアニメーションのペースに合わせて少しずつ実行していくようにするといいだろう．

例をみてみよう．

```python
import random
import simpy
from p5 import *

def manager(env):
    env.model.ordered = False  # no back order to receive
    env.stocktake = env.event()  # create the first signal (event)
    while True:
        yield env.stocktake
        report(env)
        if not env.model.ordered and env.model.level <= 10:
        # only when no back order to receive
        # reorder point = 10
            env.process(deliverer(env))  # activate deliverer
            env.model.ordered = True  # back order will be received
        env.stocktake = env.event()  # create the next signal (event)

def deliverer(env):
    yield env.timeout(5)  # delivery lead time = 5
    env.model.put(20)  # back order is recieved
    env.model.ordered = False  # no back order to receive

def customer(env):
    while True:
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        how_many = random.randint(1, 3)
        env.model.get(how_many)
        env.stocktake.succeed()  # signal for stocktaking (event)

def report(env):
    print('[{}] current level: {}, orderd: {}, queue length: {} '.format(round(env.now), env.model.level, env.model.ordered, len(env.model.get_queue)))

def simpy_setup():  # setup of simpy environment
    env = simpy.Environment()
    env.model = simpy.Container(env, init=10)  # model is marely a Container
    env.process(manager(env))
    env.process(customer(env))
#    env.run(until=200)
    return env

# ---------- code for visualization ----------
def setup():  # setup of p5py sketch
    size(600, 600)

def draw():  # drawing loop of visualization by p5py
    background(255)
    if frame_count <= 200:
        simpy_env.run(until=frame_count)  # keep pace with draw()
        level = simpy_env.model.level
        queue = len(simpy_env.model.get_queue)
        ordered = simpy_env.model.ordered
    else:
        exit()

    translate(250, 350)
    line((-200, 0), (300, 0))
    if queue > 0:
        fill(127, 0, 0)
        rect((0,0), 100, 10 *queue)  # how many customers are waiting
    if level > 0:
        fill(127)
        rect((0,0), 100, -10 *level)  # how many items at hand
    translate(0, -10 *level)
    if ordered:
        fill(255)
        rect((0,0), 100, -10 *20)  # size of back order to arrive
# ---------- ---------- ---------- ----------

if __name__ == "__main__":
    simpy_env = simpy_setup()
    run(frame_rate=10)  # you can control the speed here
```

もとの例にあった`main()`関数を微修正した，`simpy_setup()`という名称の関数を用意した．これは，SimPy側の設定を担当する関数であり，環境`env`を返す．この戻り値をsimpy_envという変数で受け取って，p5py側で利用している．

アニメーションの内容は実際に実行して確認してみてもらえたらと思うが，在庫量とバックオーダーが積み上げ型の棒グラフで表示され，欠品が出たときは，待機中の顧客の人数が下方向の棒グラフでわかるようになっている．

# まとめ

今回はSimPyで実行しているシミュレーションの途中経過をアニメーション表示する方法の例を紹介した．ここで紹介した以外にも様々な方法が考えられると思う．コンソールに文字が表示されるだけでは味気ないので，ぜひアニメーションにも挑戦してみてほしい．次回はシミュレーションにエージェントを追加する例を簡単に紹介したい．

* [第5回へのリンク](https://qiita.com/j54854/items/644c5355a38589e4a244)
* [この記事のプログラムのソース（GitHub）](https://github.com/j54854/mySimPy)