Zenn
➿

[Package] Agents.jlを利用したマルチエージェントシミュレーション | 平面と感染症を題材として

2021/06/18に公開
Julia
#
mas
tech

Agents.jl で軽く遊んでみると面白そうだったので，パッケージのドキュメントにある「Continuous space social distancing for COVID-19」という例題を動かしながら確認していこうと思います．Agents.jlですかその設計哲学として

Simple to learn and use, yet extendable, focusing on fast and scalable model creation and evolution.

ということを目指しているらしいので (ドキュメント曰く)，気軽にやっていけると思います．

環境

おそらく必要だと思われるパッケージ群と利用したversionです (描画の裏側で他にも何かいるかもしれないですが何か適宜インストールします．上のTweetと同じ内容をやるには，OpenStreetMap関係のパッケージなども必要です)．Julia自体はJulia 1.6.0を利用しています．

  • Agents.jl v4.3.1
  • CairoMakie v0.5.4
  • Compose v0.9.2
  • InteractiveDynamics v0.15.1

Agents.jlを用いたシミュレーション

本題とは関係ないですが，Agents.jlでマルチエージェントシミュレーションを実装するときの構造について最初にまとめておきます．

概要

一般的なマルチエージェントシミュレーション (MAS) については世の中の本や論文などを読んでください．大まかには以下の要素を持っている系と，その系を時間発展させる仕組みを用いて MAS と呼ばれるものが実現されていると考えていいと思います．

  • 存在するもの
    1. エージェントと呼ばれる主体
    2. エージェントなどが存在する環境
  • 動作の仕組み
    1. 周りの環境と情報などをやり取りしながらエージェントが行動する
    2. エージェントの行動が環境と他のエージェントに影響を与える
    3. ある離散時間ステップで行動や観測が行われる

例えば道路環境の MAS を考えると，エージェントとしては歩行者や車両 (それぞれ目的地を持っている)，環境としては道路などのエージェントたちを囲んでいる世界，エージェントの行動は移動に関係するもの (前進，左折，右折，停止など) に相当し，MAS を実行する主体 (私たち) は神の目の視点から定められたルールに従って行動するエージェントたちの振る舞いを観察しています．我々の目線ではエージェントの行動を動画を見るように観測することができ，時間を先に進めたり，場合によっては戻したりして，系の振る舞いを確認します．

構造

上記のような MAS を構成する要素は，Agents.jl では以下に相当します．

  • 存在するもの
    1. エージェント: Agents型のデータ (エージェントのIDや環境での位置，その他シミュレーション固有の情報を持つ)
    2. 環境: AgentBasedModel型のデータ (universal structureという名前で呼ばれる)
      • Agents.jl の様々なメソッドは，2つ (モデル(2.)とモデルに存在するエージェント(1.)) を操作するメソッドになっている
      • 例えば離散環境，グラフ環境，連続平面の環境，OpenStreetMap上の環境が実装されている
  • 動作の仕組み
    1. エージェントの操作: move_agent! や walk!，その他関数定義を行ってエージェントの振る舞いを定義する
    2. エージェントが今の環境 (モデル) において他のエージェントとどのように関係するか? として，周辺のエージェントのIDを取得したり，それらのデータ (struct data) を変更したりする．interacting_paris や nearest_neighbor，nearby_ids などモデルに応じて様々な方法で「周辺のエージェント」を取得できるようになっている．
    3. 時間ステップは step! や agent_step!，model_step! のような名前のメソッド (中身が変化するので!がつきます) を実装する．これらを用いて Agents.jl が時間発展と描画を担当する．

Agents.jlと感染症の例題 (元題: Continuous space social distancing for COVID-19) | Agents.jl を利用した MAS入門パート

ここからはドキュメントに従って，MAS の構成要素であるエージェントや系，時間発展について記述していき，動作させてみることにしましょう．イメージとしては次のようなシナリオです．

  • エージェントは2次元平面状を，各エージェントの定めた方向に定めた速度で直線で移動する
  • ビリヤードのように衝突した場合に反転する
  • SIRモデルに基づいて感染が発生する

ちなみにこの MAS についてはこちらの記事; Why outbreaks like coronavirus spread exponentially, and how to “flatten the curve” から部分的に影響を受けている (partially inspired) とのことです．

エージェント

エージェントは AbstractAgent 型を引き継ぎ，idと位置posを持っています．今二次元平面を想定しているので位置は2次元平面の位置で，ついでに移動速度と重さ (mass) を設定しておきます．ちなみにサンプルなどを確認したところ，エージェントの定義はマクロを使って実施することが多いみたいです．

using Agents, Random

mutable struct Agent <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
end

モデル

2次元平面を移動させるので，ContinousSpaceを利用します．連続のように見えているのですが適当なステップ幅で管理されているっぽいです．今periodic = falseがデフォルトで設定されているので，端が循環しているような計算がされます (Agents.jlが勝手にやります)．領域としては [0, 1] \times [0, 1] の平面です．

モデルに対して上で作成したエージェントを 500体 追加します．ここで位置は適当な場所 (ランダム)，速度はランダム角度に対して sincos を計算して作った適当なベクトルに定数倍 (speed) したもので，とりあえずmassは1.0とします．add_agent!の引数順序が少し気持ち悪いのですが，これは関数呼び出しとしてはoptionalな形です (例題だとなぜかこの形が使われていました)．こちらにいろんな方法でのagent追加があります．

function ball_model(; speed=0.002, n_agents=500)
    space2d = ContinuousSpace((1, 1), 0.02)
    model = ABM(Agent, space2d, properties = Dict(:dt => 1.0),
                rng = MersenneTwister(42))

    for ind in 1:n_agents
        pos = Tuple(rand(model.rng, 2))
        vel = sincos(2π * rand(model.rng)) .* speed
        add_agent!(pos, model, vel, 1.0)
    end
    return model
end

例1. 平面を移動するエージェント

これでエージェントと，エージェントを追加したモデルを用意したので，時間発展をさせて描画してみます．エージェントの時間発展として agent_step! を実装しますが，ここでは既に用意されている move_agent! というメソッドを使うことにします．名前の通りエージェントを与えられたモデルの上で，モデルの設定した時間幅 dt 分だけ時間発展してくれます．

描画には abm_video を利用し，描画するためにモデルと時間発展のメソッドを渡します．後ろの方は描画のパラメータで，こちらに記述がありますが描画のスピードに対するものです．

agent_step!(agent, model) = move_agent!(agent, model, model.dt)

using InteractiveDynamics
using CairoMakie

abm_video("output1.gif", model, agent_step!;
          frames=50, spf=2, framerate=25, resolution=(200, 200))

Zennのgifファイルが1.5MB以下のため小さくしてあります (200x200)．

ビリヤードっぽいエージェントのインタラクション

インタラクションですが，自分で実装しても良い (はず) ですが，ここでは Agents.jl の内部で実装されている elastic_collision! を実装します．弾性衝突のイメージはこちらのページを参照してください．あるモデルにおいて，十分に近い2体のエージェントは相互に影響を及ぼします．十分に近いエージェントを判定するためには， interacting_pairs を用いて計算します．この計算とエージェントの質量 (mass) を用いてモデルを更新します．モデルを更新するための関数として model_step! という名前で以下を定義します．

function model_step!(model)
    # 十分に近い2つのエージェントは弾性衝突する
    for (a1, a2) in interacting_pairs(model, 0.012, :nearest) 
        elastic_collision!(a1, a2, :mass)
    end
end

エージェントの移動とモデル内のエージェント間の相互作用を用いて動画を作成します．単純にagent_step!だけじゃなくmodel_step!を渡します．

model = ball_model(n_agents=100)
abm_video("output2.gif", model, agent_step!, model_step!;
          frames=50, spf=1, framerate=25, resolution=(400, 400))

さすがに200x200だと分からなかったのでエージェント数を減らし，スピードを遅くして描画してみました．目を凝らしてみると，当たったエージェント同士の起動がクッと変わっていることが分かりますね．Agents.jl自体はビリヤード計算用のソフトではないので，この弾性衝突計算は完全なものではないとドキュメントに書いてありますので注意しましょう．

全員動くのが面倒なので動かないエージェントを追加する

質量が無限で速度が0のエージェントを追加すれば良いです．今，Agentの型はmutableになっているので，500のうち400を固定するために質量と速度を変えるには次のように実装すれば良いです．

model = ball_model()
for id in 1:400
    agent = model[id] # エージェントはIDを持っています
    agent.mass = Inf
    agent.vel = (0.0, 0.0)
end

このモデルを同じように描画してみましょう．500だと多すぎたので300のうち200を無限質量にしてみました．

Agents.jlと感染症の例題 (元題: Continuous space social distancing for COVID-19) | 感染症パート

ここからはSIRモデルを思い浮かべつつ，実装を見ていきます．

SIRを考慮するエージェント

基本的なデータはAgentと同じですが，感染してからの日数と状態(S or I or RをJuliaのSymbol :S, :I, :Rで持ちます)，またついでにSからIへ変遷する確率βを持ちます．またIからRへの変遷については，death_rate というパラメータを考えて，if文で変遷することにします (後ほど)．

まずはエージェントのデータ構造です．

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
    days_infected::Int  # number of days since is infected
    status::Symbol  # :S, :I or :R
    β::Float64
end

SIRを考慮するモデル

まず次のような病気の前提を頭で思い浮かべます．

  • 1日はシミュレーションステップ24ステップとする (1ステップ1時間)
  • 1度感染すると30日間感染する (infection_period)
  • 感染が発覚するまでに14日間必要 (detection_time)
  • 5%の確率で再感染する (reinfection_probability)
  • 4.4%の確率で亡くなる (death_rate)
  • 初期に1000人中5人が感染していた
  • S => Iへの遷移確率は [0.4, 0.8] の間とする

上が実装されたコードが以下のものです (調整済み)．

function sir_initiation(;
    infection_period = 30 * steps_per_day,
    detection_time = 14 * steps_per_day,
    reinfection_probability = 0.05,
    isolated = 0.0, # in percentage
    interaction_radius = 0.012, dt = 1.0, speed = 0.002,
    death_rate = 0.044, # from website of WHO
    N = 1000, initial_infected = 5, seed = 42, βmin = 0.4, βmax = 0.8,
)

    # パラメータと領域
    properties = Dict(
        :infection_period => infection_period,
        :reinfection_probability => reinfection_probability,
        :detection_time => detection_time,
        :death_rate => death_rate,
        :interaction_radius => interaction_radius,
        :dt => dt
    )
    space = ContinuousSpace((1,1), 0.02)
    model = ABM(PoorSoul, space, properties=properties, rng=MersenneTwister(seed))

    # Nエージェントの追加
    for ind in 1:N
        pos = Tuple(rand(model.rng, 2))
        status = ind ≤ N - initial_infected ? :S : :I
        isisolated = ind ≤ isolated * N
        mass = isisolated ? Inf : 1.0
        vel = isisolated ? (0.0, 0.0) : sincos(2π * rand(model.rng)) .* speed
        β = (βmax - βmin) * rand(model.rng) + βmin
        add_agent!(pos, model, vel, mass, 0, status, β)
    end

    return model
end

このモデルを描画だけしたものが公式サイトにあります．

SIRを考慮するモデルの時間発展

エージェントとモデルを作成したので，時間発展を定義します．基本的な考え方 (何が実装されているか) については

  • エージェントのペア a1, a2 について，SIR的な考え方で状態 (:S, :I, :R) が変化する (モデルのパラメータ次第)
  • モデルの変遷では全エージェントが弾性衝突するように移動する
  • エージェント単体の変化として，病気からの復帰・死亡があります

です．

S→Iの遷移

再感染の確率がある場合，それを考慮して状態遷移 (:S, :I, :R) します．

function transmit!(a1, a2, rp, model)
    # for transmission, only 1 can have the disease (otherwise nothing happens)
    count(a.status == :I for a in (a1, a2)) ≠ 1 && return
    infected, healthy = a1.status == :I ? (a1, a2) : (a2, a1)

    rand(model.rng) > infected.β && return

    if healthy.status == :R
        rand(model.rng) > rp && return
    end
    healthy.status = :I
end

弾性衝突部分 (+ S→I)

上で実装したものと同じですが，弾性衝突する際にS→Iの判定を行い，既に感染している人から他の人へうつる可能性があります．

function sir_model_step!(model)
    r = model.interaction_radius
    for (a1, a2) in interacting_pairs(model, r, :nearest)
        transmit!(a1, a2, model.reinfection_probability, model) # ここがS→I
        elastic_collision!(a1, a2, :mass)
    end
end

エージェント単体の変化

実装の中身は移動することと状態を更新すること (感染からの日数をカウントすること)，そして一定の確率でお亡くなりになる (modelから消えるので kill_agent) ことです．治った場合はリセットです (0日目)．

function sir_agent_step!(agent, model)
    move_agent!(agent, model, model.dt)
    update_agent!(agent)
    recover_or_die!(agent, model)
end

update_agent!(agent) = agent.status == :I && (agent.days_infected += 1)

function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand(model.rng) ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end

描画した結果

これまでコピペしながら実装してきた内容を描画したものがこちらになります．

公式のコピペだと微妙に動かない (transmit!) ので model を追加で渡して雰囲気で直してください．いい感じに感染が広がって赤くなっていくのが分かりますね (このステップ数だと，見た目だけでは治っていくところまで見えない気がします)．こちらの例はN=300にした上で，GIFの容量を抑えるためにギュッと小さくして描画したものなので，大きいものを見たい場合には手元で動かしてみてください．

統計量の観察

ここまでで見た目が面白いgifファイルが作れたので，感染者数などの量を計測したいと思います．Agents.jl が測定してDataFramesで返す仕組みを持っているので簡単です．

以下に感染者数と治った人の数を測定するメソッドと，測定したいよ～～～という気持ち (Agentの:statusに対してinfectedやrecoveredを適用したいんだよという気持ち) を用意します．

infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)
adata = [(:status, infected), (:status, recovered)]

測定するには，今まで描画 abm_video していたものを run! というメソッドを用います．引数に上の「測定したいよ～～～」という気持ちのタプルと，描画するときと同じようにモデルと時間遷移するためのメソッドを突っ込みます．以下のコードはいろんなパラメータを変えながら測定するコードは次のようなものです (コピペ)．

r1, r2 = 0.04, 0.33
β1, β2 = 0.5, 0.1
sir_model1 = sir_initiation(reinfection_probability = r1, βmin = β1) # 設定1
sir_model2 = sir_initiation(reinfection_probability = r2, βmin = β1) # 設定2
sir_model3 = sir_initiation(reinfection_probability = r1, βmin = β2) # 設定3

data1, _ = run!(sir_model1, sir_agent_step!, sir_model_step!, 2000; adata) # 測定1
data2, _ = run!(sir_model2, sir_agent_step!, sir_model_step!, 2000; adata) # 測定2
data3, _ = run!(sir_model3, sir_agent_step!, sir_model_step!, 2000; adata) # 測定3

# 末尾を見る
data1[(end-10):end, :]

末尾を見たデータフレームの中身を見ると，世代のステップ数と各統計量が取れていることが分かりました (例えば infected_status が上で指定していたエージェントのstatusにinfectedを適用して測定した結果に相当します)

CairoMakie.jl を利用してグラフにするとこうなります．

(おまけ) ソーシャルディスタンスの効果を見積もる

これまでにSIRモデルを考慮したマルチエージェントシミュレーションのモデルを動かしながらAgents.jlについて雰囲気を学んできました．

ところで，エージェントを作ったときに動かないエージェント (質量無限，移動0) を実装していたことを思い出しましょう．今回はSIRを考慮したモデルと動かないエージェントを用いて，「ソーシャルディスタンス」（最近聞かなくなりましたが）についてシミュレーションすることにします．

全体のうち 80\% の人が引きこもって移動しなくなった謎の二次元平面を考えシミュレーションして計測を行います．

r4 = 0.04
sir_model4 = sir_initiation(reinfection_probability = r4, βmin = β1, isolated = 0.8)

# あとはCairoMakie.jlでグラフを重ねる

行った計測 (4番目) を，これまでに作成した3つの計測のグラフに重ねてみます．

面白い結果になりましたね (?)．最後にAgents.jlのドキュメントの一文を引用して終わります．

animations are cool, but science is even cooler.

Discussion