Dynamic Models of Segregation (Schelling 1971)をリアルタイムに表現するコードを更新したのでそのメモを残しておきます．
今回はNetLogo風に，エージェントの幸福度をリアルタイムに計算してプロットする部分を追加してみました．
Mathematica ノートブックは以下のサイト内で公開しているseg.nbというファイルをDLしてください．
使い方は，ノートブック内にも書いていますが，
- Module セルを読み込む
- 評価 [RightArrow] 動的更新の実行にチェックを入れる
- Manipulate セルを読み込む
- resetボタンを押す
- run simulationのチェックボックスをon
- 終了すると END!! が出て自動的に停止する
- 条件を変えて4.から繰り返し
です．最終的なアウトプットはこんな感じです．
- 赤と青の個体がトーラス上をランダムに移動します
- 各個体は自分の近傍内に同色の仲間が35%以上存在すると，その位置に定着します
- 左のバーは平均的な幸福度（周りの同色エージェント割合）です
- このシミュレーションは，強い排外意識がなくても，分居が自然に生じる様子を表しています
ノートブックを読み込む段階で動的更新をonにしていると，このUIが表示されるはずですが，エラーが出る場合もあるので 最初はオフにしておき，上記の手順で動的更新を許可すれば問題ないはずです．
このnbはWolfram Demonstrations Projectで公開されているコードを参考にして作りました．
特に，GarbageCollectionByAntsというコードを参考にさせていただきました．
使用関数の定義
まずエージェントの初期配置を決める関数．size×sizeの行列上の番地でエージェントの位置を指定する． 赤エージェントの数は，rednで，赤エージェントの数は，bluenで指定する．
(*赤agent 青agentの初期配置*) initialset[size_, redn_, bluen_] := Module[{a, b, n = redn + bluen}, Developer`ToPackedArray[(* Integer型指定で高速化 *) a = RandomSample[Range[size*size], n];(* 総エージェント数だけ1次元配列から非復元抽出 *) b = Flatten[Table[{i, j}, {i, 1, size}, {j, 1, size}], 1];(* 2次元座標 *) (* 抽出した1次元配列番号を座標に変換する． 前半は赤，後半は青用 *) {Table[b[[a[[i]]]], {i, 1, redn}], Table[b[[a[[i]]]], {i, redn + 1, n}]}, Integer] ]; (* 非復元抽出で初期配置を与える． output: {{赤の位置情報}，{青の位置情報}} の順で出力する initialset[size,redn,bluen][[1]] が 赤初期配置 initialset[size,redn,bluen][[2]] が 青初期配置 *)
初期グリッドに0を並べる手続きです．動的表現を使わない場合は，わざわざ定義する必要はありません．
(*初期グリッド*) zerogrid[size_] := Table[0, {size}, {size}]
エージェントの位置情報を色に変換する関数です．この関数の引数である位置情報がDynamicで更新されるのでエージェントが動いているように見える，という仕組みです．
(* agentの位置情報を検索して，可視化する関数 *) visualize[{redposi_, blueposi_}, size_] := Module[{zero1 = zerogrid[size](*Mnipulate内で指定*)}, (* Embedding agents into zero1 for visualization 引数のredposi_,blueposi_はそれぞれ「赤」と「青」の現在位置 *) Do[zero1[[ Sequence @@ redposi[[i, {1, 2}]]]] = 1, {i, Length[redposi]}]; Do[zero1[[ Sequence @@ blueposi[[i, {1, 2}]]]] = 2, {i, Length[blueposi]}]; ArrayPlot[zero1, ColorFunction -> (If[# == 0, White, If[# == 1, Darker[Red], Darker[Blue]]] &), ColorFunctionScaling -> False, ImageSize -> {300, 300}, AspectRatio -> Automatic, Frame -> False, Mesh -> All, Axes -> False] ];
全セル（トーラス空間）のムーア近傍の座標を出力する関数です．エージェントの位置毎に計算するよりも，グローバルに一回計算して，必要があるやつだけ呼び出したほうが早いんじゃないかと思って定義した関数です．計算は一回きりなので，メモリの節約にはなっていると思うのですが，実際に比較はしてません．ただ，コードの構造上はこの方法が理解しやすいと思います．
(* トーラス上に配置されたエゴのムーア近傍を出力する関数 *) moore[t_(* agent position *), size_] := Module[{m = size, p1, p2, p3, p6, p4, p7, p8, p9, network}, p6 = Table[ If[t[[i]][[2]] != m, t[[i]] + {0, 1}, t[[i]] + {0, -(m - 1)}], {i, Length[t]}]; p9 = Table[Which[(* +-している1を一般化すれば近傍の距離を自由に定義できる*) (t[[i]][[1]] != 1) && (t[[i]][[2]] != m), t[[i]] + {-1, 1}, (t[[i]][[1]] == 1) && (t[[i]][[2]] == m), t[[i]] + {m - 1, -(m - 1)}, (t[[i]][[1]] == 1) && (t[[i]][[2]] != m), t[[i]] + {m - 1, 1}, (t[[i]][[1]] != 1) && (t[[i]][[2]] == m), t[[i]] + {-1, -(m - 1)}], {i, Length[t]}]; p7 = Table[Which[ (t[[i]][[1]] != 1) && (t[[i]][[2]] != 1), t[[i]] + {-1, -1}, (t[[i]][[1]] == 1) && (t[[i]][[2]] == 1), t[[i]] + {m - 1, (m - 1)}, (t[[i]][[1]] == 1) && (t[[i]][[2]] != 1), t[[i]] + {m - 1, -1}, (t[[i]][[1]] != 1) && (t[[i]][[2]] == 1), t[[i]] + {-1, (m - 1)}], {i, Length[t]}]; p3 = Table[Which[ (t[[i]][[1]] != m) && (t[[i]][[2]] != m), t[[i]] + {1, 1}, (t[[i]][[1]] == m) && (t[[i]][[2]] == m), t[[i]] + {-(m - 1), -(m - 1)}, (t[[i]][[1]] == m) && (t[[i]][[2]] != m), t[[i]] + {-(m - 1), 1}, (t[[i]][[1]] != m) && (t[[i]][[2]] == m), t[[i]] + {1, -(m - 1)}], {i, Length[t]}]; p1 = Table[Which[ (t[[i]][[1]] != m) && (t[[i]][[2]] != 1), t[[i]] + {1, -1}, (t[[i]][[1]] == m) && (t[[i]][[2]] == 1), t[[i]] + {-(m - 1), (m - 1)}, (t[[i]][[1]] == m) && (t[[i]][[2]] != 1), t[[i]] + {-(m - 1), -1}, (t[[i]][[1]] != m) && (t[[i]][[2]] == 1), t[[i]] + {1, (m - 1)}], {i, Length[t]}]; p4 = Table[ If[t[[i]][[2]] != 1, t[[i]] + {0, -1}, t[[i]] + {0, (m - 1)}], {i, Length[t]}]; p2 = Table[ If[t[[i]][[1]] != m, t[[i]] + {1, 0}, t[[i]] + {-(m - 1), 0}], {i, Length[t]}]; p8 = Table[ If[t[[i]][[1]] != 1, t[[i]] + {-1, 0}, t[[i]] + {(m - 1), 0}], {i, Length[t]}]; network = Developer`ToPackedArray[ Table[{p6[[i]], p4[[i]], p2[[i]], p8[[i]], p9[[i]], p7[[i]], p3[[i]], p1[[i]]}, {i, Length[t]}], Integer] ];
エージェントの幸福度を計算する関数です．最小値0，最大値1です．位置情報と幸福度をセットにして出力します．
(* ムーア近傍情報に基づき幸福度（ムーア近傍内の同色エージェント割合）を計算, 位置情報に幸福度を追加 引数 agentposiは単色エージェント *) happy[agentposi_, size_] := Module[{neloca, happy3}, neloca = moore[agentposi, size]; happy3 = Table[Plus @@ Table[If[MemberQ[agentposi, neloca[[i]][[j]]], 1, 0], {j, 8}], {i, 1, Length[neloca]}]; happy3 = happy3/8;(* 同色エージェント数を *) Table[Append[agentposi[[i]], happy3[[i]]], {i, Length[agentposi]}] ]; (*入力は {{1,1},{2,1},{3,1}} のような単色エージェントの位置情報．出力は {{1,1,0.2},{2,1,0.3},{3,1,0.8}} のような 位置情報+幸福度 *)
最後にエージェントの移動を管理する関数です．設定した同色エージェント割合（閾値）に達していない場合に移動します．
(* 幸福情報に基づき，移動する visualizeの引数を生成する *) move[agents_, size_, redth_, blueth_, movern_] := Module[{preloca, zerolocate, red0, blue0, red1, blue1, k1, k2, AllocationforRed, l1, l2, AllocationforBlue}, red0 = agents[[1]]; blue0 = agents[[2]]; red1 = happy[red0, size]; blue1 = happy[blue0, size]; preloca = Flatten[{red0, blue0}, 1];(* 埋まっているセル *) zerolocate = Complement[Tuples[Range[size], 2], preloca];(*ゼロセルの番地リスト*) (*red agentからmovern人を選び動かす*非復元抽出* *) k1 = RandomSample[Range[Length[red0]], movern]; k2 = Select[k1, red1[[#, 3]] < redth &]; AllocationforRed = RandomSample[zerolocate, Length[k2]]; red[[k2]] = AllocationforRed; (*Blue agentからmovern人を選び動かす*非復元抽出* *) l1 = RandomSample[Range[Length[blue0]], movern]; l2 = Select[l1, blue1[[#, 3]] < blueth &]; AllocationforBlue = RandomSample[Complement[zerolocate, AllocationforRed], Length[l2]]; blue[[l2]] = AllocationforBlue; {red, blue} ];
さて，以上の関数をforで回すのではなく，Manipulateの中でRefreshします．
するとリアルタイムでエージェントがピコピコ動きます．このような動的表現が，MathematicaのDynamic関数によって実行できます．
以下Manipulateの中身です．
(* ver110:幸福度をlistplotでリアルタイムに更新 自動停止機能を追加，happy level モニタリングをDynamic追加 *) Manipulate[ Refresh[ If[moving, {red, blue} = move[{red, blue}, s, redth, blueth, nomove]](* 更新している部分 *); If[(Min[happy[red, s][[All, 3]]] >= redth) && (Min[happy[blue, s][[All, 3]]] >= blueth), moving = False; state = Text[Style["END !", Red, 20]], moving = moving]; state2 = Mean[happy[red, s][[All, 3]]] // N; state3 = Mean[happy[blue, s][[All, 3]]] // N; (* agentのhappy を 監視して，全員が満足したらmoving をFalseに切り替え, Dynamic[state]表示をendに変える． moving の値がFalseになると {red,blue}の更新が終わり， "run simulation"チェックボックスがoffになる*) visualize[{red, blue}, s], UpdateInterval -> If[moving, 0.01, Infinity]](* reflesh end *), {{s, 20, "length of a side"}, 10, 30, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{redn, Floor[2*s], "number of red"}, 1, Floor[(s*s)/2], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{bluen, Floor[2*s], "number of blue"}, 1, Floor[(s*s)/2], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{redth, 0.35, "red threshold"}, 0, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{blueth, 0.35, "blue threshold"}, 0, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{nomove, Floor[(s)/5], "number of mover"}, 1, Min[{redn, bluen}], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{moving, False, "run simulation"}, {True, False}}, Button["reset", moving = False; state = "going..."; x = initialset[s, redn, bluen]; red = x[[1]]; blue = x[[2]], ImageSize -> Medium], Dynamic[{"state=", state}], Dynamic[ListPlot[{{{1, state2}}, {{2, state3}}}, Filling -> Axis, PlotRange -> {{0, 3}, {0, 1}}, PlotMarkers -> {Red, Blue}]], SaveDefinitions -> True, TrackedSymbols :> {moving, red, blue, redn, bluen, s}, SynchronousUpdating -> True, ControlPlacement -> Left ]
なにゆえ，このManipulateでエージェントが動くのかを，ざっと解説します．
Manipulateの直接の実行対象となっている関数は Refresh[ ]です．Refresh[ ]内で，上で定義した一連の関数をループさせています．
ループさせている関数の中には，visualizeがあるので，Refreshされる度に，グラフィック出力が変化します．このRefreshするタイミングを
UpdateInterval ->で指定します．intervalが短い（上記コードでは0.01）ので，CPUさえ高性能なら，ヌルヌル動くというわけです*1．
実験としてスライダの下部にDynamic objectを2つ埋め込みました．一つはシミュレーションが停止したかどうかを示すstateです．
停止条件を満たすとここがEND！と表示されます．
もう一つが，今回修正したエージェントの幸福度（リアルタイム）です．色別の幸福度の平均値を計算してListPlotする命令をDynamicに放り込むだけでうまくいきました．
Dynamicは，誤って使用するとPCに負担がかかりすぎてカーネルが落ちるため，なかなか使いこなすのが難しい関数ですが，うまくいけば動的な表現の幅が広がるので是非
使ってみてください．
（今みなおすと，まだ無駄な計算が残っており，もっと効率化できそうです．どのあたりが修正できそうか，気づいた学生さんは授業時間内に教えてください）
*1:厳密に言えば，このオブジェクトは最速で更新されるために，インターバルの指定は無視されます．遅延評価する実行文の場合，インターバルは指定通りになります．これはRefreshの仕様です
