Wolfram's Cellular Automaton(CA)

code: https://colab.research.google.com/drive/1-W0CJG2HYR8X5dWuBigiW8hRmH3mIVI7

参考：Mathematica複雑系のシミュレーション14章 「複雑性」

Steven Wolframの一次元CA

周期的境界条件を持つ1行の格子サイトに対し、確率的で均一（一様）な局所ルールを適用した場合、発展形態として次の4種類に落ち着くことをウルフラムは発見しました。

  • 均一状態：集積点limit pointに似る
  • 周期的構造：極限軌道 limit cycleに似る
  • カオス的パターン：strange attracterに似る
  • 複雑な局所構造

この場合、最後の「複雑な局所構造」はそれ以外の3つと異なり、物理モデルの定式化が難しく、コンピュータ・シミュレーションを必要とします。逆にいえば、定式化された物理モデルを持たない数理科学現象の時間発展を簡単な準備だけでシミュレーションするのにCAは向いているといえます。一方、均一状態、周期的構造、カオス的パターンなどの物理現象の可視化を容易に行える力をCAは持っています。

最初に、最も容易で、かつその時間変化が自明ではないCAの例を一つ挙げます。

In [ ]:
{{1} -> 0, {0} -> 0}
{{1} -> 0, {0} -> 1}
{{1} -> 1, {0} -> 0}
{{1} -> 1, {0} -> 1}

サイトの値（カラーk）は0か1の2種類（k=2）、近傍の値は考慮せず自身の値のみを想定（r=0）し、たとえば最初のルール{{1}->0,{0}->0}が適用された場合、値1を持つサイトは0を、値0を持つサイトは0となります。

次に k=2, r=1のとした時、各サイトは両側に値0か1を持つ最近接サイトが来ることになるので、この系の可能な近傍の値は次のようになります（中央のサイト値1の組み合わせ4通り、中央のサイト値0の組み合わせ4通り）。

In [ ]:
{{1, 1, 1}, {1, 1, 0}, {1, 0, 1}, {1, 0, 0},
 {0, 1, 1}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}

この系に関し、ルールを一つ考えます。

In [ ]:
{{1, 1, 1} -> 1, {1, 1, 0} -> 1, {1, 0, 1} -> 1, {1, 0, 0} -> 0,
 {0, 1, 1} -> 1, {0, 1, 0} -> 0, {0, 0, 1} -> 0, {0, 0, 0} -> 1}

では、こうしたルールの数は何種類になるでしょう。各組み合わせ（$8通り=2^3$）に対し、1か0を取るとすれば、ルールセットの総数は$2^8=256$個となります。したがって一般的なルールセット総数は

In [ ]:
totalRuleSet[k_, r_] := k^(k^(2 r + 1))

となります。上記例でいえば、

In [ ]:
totalRuleSet[2, 0]
Out[ ]:
4
In [ ]:
totalRuleSet[2, 1]
Out[ ]:
256

と計算できます。たとえば k=3（各サイトが3種類の値を取る）、r=1（両側の近接サイト数=1）の場合、想定ルール数は次のような巨大な数になります。

In [ ]:
totalRuleSet[3, 1]
Out[ ]:
7625597484987

次にk, r, w（ルール番号）からルールセットを生成する関数 ruleSet を定義してみます。最初にルール番号 w について確認しておきます。たとえば、上記ルールセット例

{{1, 1, 1} -> 1, {1, 1, 0} -> 1, {1, 0, 1} -> 1, {1, 0, 0} -> 0, {0, 1, 1} -> 1, {0, 1, 0} -> 0, {0, 0, 1} -> 0, {0, 0, 0} -> 1}

を例に取ると、右辺だけを並べると 11101001 となるので、これを k (=2) を基準に十進数に直し、それをルール番号 w とします。

In [ ]:
w = 2^^11101001
Out[ ]:
233

逆を求めるには IntegerDigit関数を利用します。

In [ ]:
IntegerDigits[w, 2]
Out[ ]:
{1, 1, 1, 0, 1, 0, 0, 1}

ruleSet

In [ ]:
ruleSet[k_, r_, w_] := Module[{fillCode, tuplets},
  fillCode[x_, y_, z_] := 
     Join[Table[0, {x^(2 y + 1) - Length[#]}],
        #] &[IntegerDigits[z, x]];
  tuplets[p_, s_] :=
     Reverse[Flatten[Apply[Outer,
        Prepend[Table[Range[0, p - 1], {(2 s + 1)}], List]], 2 s]];
  	MapThread[Rule, {tuplets[k, r], fillCode[k, r, w]}]
  ]

ルール番号0,1,2,３それぞれの場合のルールセットを生成します。

In [ ]:
Table[With[{k = 2, r = 0},
   ruleSet[k, r, w]], {w, 0, 3}] // TableForm
Out[ ]:
Output
In [ ]:
With[{k = 2, r = 1, w = 233},
 ruleSet[k, r, w]] //TableForm
Out[ ]:
Output

WolframCA

In [ ]:
Options[WolframCA] =
  {CASize -> 100,
   CATimeSteps -> 20,
   InitialState -> "Disordered"};
   
WolframCA[w_, k_, r_, opts___Rule] :=
 
 Module[{ruleSet, row, size, steps, state},
  size = CASize /. {opts} /. Options[WolframCA];
  steps = CATimeSteps /. {opts} /. Options[WolframCA];
  state = InitialState /. {opts} /. Options[WolframCA];
  row = If[state === "SingleSeed",
    ReplacePart[
     Table[0, {size}], k - 1, Ceiling[size/2]], 
    Table[Random[Integer, {0, k - 1}], {size}]];
  ruleSet = MapThread[Rule, {
     Reverse[Flatten[Apply[Outer,
        Prepend[Table[Range[0, k - 1], {(2 r + 1)}], List]], 2 r]], 
     Join[Table[0, {k^(2 r + 1) - Length[#]}], #] &[
      IntegerDigits[w, k]]}
    ];
  FixedPointList[
   (Partition[Join[Take[#, -r], #, Take[#, r]], (2 r + 1), 1]
      /. ruleSet) &,
   row, steps]]
In [ ]:
ShowCA[list_List] := 
 Module[{rule, max = Max[list], colorrule},
  	colorrule = Flatten[
    Map[{# -> Hue[#/(max + 1)]} &, Range[0, max]]];
  	ArrayPlot[Reverse[list] /. colorrule]
  	]

I類の挙動

初期配置からの変化が起こり、すぐに死滅し不動になった均一状態になる。（時間軸は下から上に流れます。）

In [ ]:
ShowCA[ WolframCA[136, 2, 1] ]
Out[ ]:
Output

II類の挙動

互いにつながっていない周期的なる領域からなるパターンが出現する。

In [ ]:
ShowCA[ WolframCA[4, 2, 1] ]
Out[ ]:
Output

III類の挙動

非周期的な見かけ上カオスのようなパターンが生成される。

In [ ]:
ShowCA[ WolframCA[45, 2, 1,
    CASize -> 150, CATimeSteps -> 80] ]
Out[ ]:
Output

IV類の挙動

複雑でかなり局所的な構造を生成する。

In [ ]:
ShowCA[ WolframCA[110, 2, 1,
    CASize -> 150, CATimeSteps -> 60] ]
Out[ ]:
Output

初期配置を一点に限った場合。

In [ ]:
ShowCA[ WolframCA[30, 2, 1,
    CASize -> 400, CATimeSteps -> 300, 
    InitialState -> "SingleSeed"] ]
Out[ ]:
Output

今まではWolfram CAの原理をMathematicaのコードのみで書いてきました。現在は、上記関数はMathematicnaの組み込み関数 CellularAutomaton として高速実装されています。次の例は上記例のCellularAutomatonによる再実装となります。

In [ ]:
ShowCA[CellularAutomaton[30, {{1}, 0}, 300]]
Out[ ]:
Output

同一ルールセットを使いランダムウォークを実現します。

In [ ]:
ListLinePlot[
 (* (-1)^: 単位ステップ、300:経過時間、{{0}}:出発時刻 *)
    Accumulate[(-1)^CellularAutomaton[30, {{1}, 0}, {300, {{0}}}]]]
Out[ ]:
Output

２次元セルラーオートマトン

9近傍 k=2 r=1 の例

In [12]:
showSites[sites_, opts___] :=
   ArrayPlot[sites, opts, Mesh -> True, PixelConstrained -> 4]
In [13]:
SeedRandom[3];
config = RandomInteger[1, {21, 21}];
showSites[config]
Out[13]:
Output

LifeGameの例

In [16]:
sites = CellularAutomaton[{224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 
       2}}}, {1, 1}}, config, 4];
showSites /@ sites
Out[16]:
Output

CellularAutomatonのListConvolveによる実装を紹介します。

In [ ]:
Clear[CAStep]
CAStep[rule_, a_] := 
 Map[rule[[18 - #]] &, 
  ListConvolve[{{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}, a, 2], {2}]
In [ ]:
sites = NestList[CAStep[IntegerDigits[224, 2, 18], #] &, config, 4];
showSites /@ sites
Out[ ]:
Output
In [ ]:
Grid[{Range[18], IntegerDigits[224, 2, 18]}, 
 Dividers -> {{}, {False, True}}]
Out[ ]:
Output
In [ ]: