IA5 - 元胞自動機

生命遊戲 (Game of Life, GoL)，取自維基百科

元胞自動機 (cellular automaton，複數automata，縮寫CA) 由John von Neumann在上個世紀50年代所提出，希望能用來模擬生物自我複製的情況，然後到了70年代由John Conway設計了生命遊戲 (Game of Life, GoL) 並刊上了科學美國人後才逐漸吸引其他研究者的目光。而在80年代由天才科學家，也是知名數學套件Mathematica®的開發者Stephen Wolfram利用統計力學的方式做了一系列的研究後，我們才對這類型的計算理論與結果有了系統性的理解。

一、元胞自動機的演化規則

CA通常由包含了\(N\)個元素的一維陣列所組成，每個元素不是\(1\)便是\(0\)。用比較圖像化的形式，將每個元素都想像成一個細胞，每個細胞不是黑\(1\)就是白\(0\)的兩種狀態，而整個陣列便代表在時間\(t\)時由\(N\)個細胞所構成的假想生物體的狀態，當然這個假想生物體就是由一系列\({\tt 101000100}\cdots\)來描述。
　現在這個生物體的狀態由時間\(t\)演化到\(t+1\)，那麼在\(t+1\)時這個生物體的每個細胞是處在\(1\)或\(0\)則由上一步的它自己與它左右兩個鄰居的狀態來決定，如下： \[ \begin{array}{cccccc} t & \cdots & {\tt 1} & {\tt 0} & {\tt 1} & \cdots\\ & & \searrow & \downarrow & \swarrow\\ t+1 & & & ? \end{array} \] 可以看見在\(t\)時三個相鄰的細胞狀態分別是\({\tt 101}\)，則處在\(0\)的那個細胞它在\(t+1\)時則由它自己原來的狀態\(0\)，以及鄰近那兩個處在\(1\)的細胞來決定。
　所以說任三個相鄰細胞的狀態，可以決定中間那個細胞在下一個時刻時的狀態，由於每個細胞可以處在\(0\)和\(1\)兩種情況，三個細胞變有\(2^3=8\)種可能的組合： \[ \begin{array}{ccccccccc} t & \tt 111 & \tt 110 & \tt 101 & \tt 100 & \tt 011 & \tt 010 & \tt 001 & \tt 000\\ t+1 & \tt 0 & \tt 0 & \tt 1 & \tt 1 & \tt 0 & \tt 0 & \tt 1 & \tt 0 \end{array} \] 上面這個例子，我取若三個細胞狀態都是\(1\) (最左)，則中間那個細胞在\(t+1\)時會變成\(0\)，同樣的若是\(\tt 110\)則中間在下一刻會變成\(0\)⋯⋯ ，由於二進位代碼\(\tt 00110010\)對應到十進位的\(50\)，Stephen Wolfram便將元胞自動機遵從上面的演化規則稱作Rule 50，另外用前面提到的黑白來表述的話，則如圖1所示。

圖1. Rule 50

　由於在任意時刻\(t\)任三個相鄰的細胞總共有\(8\)種可能的情況，而中間那個細胞在下一個時刻不一定要強制是Rule 50，比如說圖1倒數第二個原來是\(\tt 001\)而在下個時刻中間那個細胞我可以讓他變成\(0\)而不是Rule 50的\(1\)，那麼這便產生了新的規則\(\tt 00110000\)，而這組二進位代碼對應到48，這個新的規則就叫做Rule 48。總括來說，\(8\)種相鄰的情況，每種情況的中間那個細胞在下一刻可以選擇處在\(0\)或\(1\)兩個狀態，所以總共可以寫下\(2^8=256\)條規則，從Rule 0~255，Stephen Wolfram在上世紀80年便是對這256條規則所演化出來的CA做出詳細的研究。Rule 30和120的情況可以從圖2中看到。

圖2. 另外兩種CA的規則

1.1 元胞自動機的分類

Stephen Wolfram將Rule 0~255生成的CA分成四種類型：

1. 大部分的初始狀態都迅速演變成一穩定 (stable) 且均勻 (homogeneous) 的圖像，初始的隨機性不再
2. 大部分的初始狀態都迅速演變成一穩定且具有週期性規律震盪 (oscillating) 的圖像。初始條件的隨機性有部分被保留但部分在演化的過程中會消失，這些消失的隨機性會被CA本身的規律所取代，並形成局部 (local) 的圖像，這些局部的圖像在往後的演化仍會是局部的
3. 大部分的初始狀態都會迅速變成一偽隨機或混沌的狀態，任何有規律的圖像都會快速地被埋沒在周圍的噪聲當中，這些導致局部改變的情況會隨著演化擴散到整幅圖像當中
4. 大部分初始狀態會迅速演變成特定的結構，並與周遭的環境產生交互作用，這些結構會維持極長一段時間。據推測第4類最終會回歸到第2類的狀態，但可能需要非常長的時間來達成

以上所謂的初始條件都是給定比如初始陣列是只有中間的細胞狀態是\(1\)，其它都是\(0\)的狀況，接著我看最後的演化會變成什麼樣子，我們將在下一節的時候再回到這裡來。

二、實作元胞自動機

這裡python的作法參考自Allen Downey的Think Complexity，感謝Prof. Downey寫了這麼有趣的一本書！
　首先我們需要一條將十進位0~255的數字轉成二進位的方式，因為這個八碼的規則會決定每個細胞在下一個時刻的狀態 (當然是受到這個時刻它和它鄰居之間的狀態所影響)。這可以使用np.unpackbits來達成，該函數接受任何元素為8位元長度的無號整數 (unsigned interger) 陣列並轉成二進位代碼。所以我們可以建立規則轉換函數

# Generate Stephen Wolfram's rule in binary form
def rule_tab(num):
    num = np.array([num],dtype=np.uint8)
    return np.unpackbits(num)[::-1]

故rule_tab接受任何介於0~255之間的數num，並將它轉成8碼的二進位規則，至於最後為什麼要加[::-1]讓輸出的順序顛倒會在2.1節說明。
　接下來便是怎麼讓python去利用這條規則來計算下一步每個細胞的狀態。我們先給定一個生物體擁有\(9\)個細胞，初始狀態是\(\tt 100101101\)，然後我們想依照Rule 50 (圖1) 算到第\(20\)步時這個生物體每個細胞的狀態，我們可以先用np.zeros先建立一個\(20\times 9 \)的零矩陣，接下來將這個矩陣的第0列替換成\(\tt 100101101\)，所以會得到 \[ {\tt array}=\left[\begin{array}{ccccccccc} {\tt 1} & {\tt 0} & {\tt 0} & {\tt 1} & {\tt 0} & {\tt 1} & {\tt 1} & {\tt 0} & {\tt 1}\\ {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0} & {\tt 0}\\ \vdots & & & & \ddots & & & & \vdots\\ {\tt 0} & & \cdots & & {\tt 0} & & \cdots & & {\tt 0} \end{array}\right].\tag{1} \] 所以array的第0列代表初始狀態，第1列代表第一步的狀態以此類推，只是此時還沒計算，一直到最後一步我們都先用\(0\)來暫時代表。
　那麼我們怎麼演化到第1步呢？我們取出array的頭兩列 \[ \begin{array}{cccccccccc} t=0: & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 1\\ t=1: & a & b & c & d & e & f & g & h & i \end{array} \] 所以在\(t=1\)的時候\(a\)~\(i\)會分別填上\(0\)或\(1\)，注意對邊界的細胞而言，像\(a\)本身只有一個右邊的鄰居，那左邊的空鄰居在計算上一個狀態時就自動補\(0\) ；而最左邊的細胞\(i\)右邊的空鄰居也是一樣的做法。所以對\(b\)細胞而言，它的上一步是\(\tt 100\)，按照Rule 50而言應有\(b=1\)，其它的細胞可以按照這個方式填上狀態。
　所以我們要寫下決定細胞狀態的狀態算法，對於生物體第\(t\)步的狀態，可以透過\(t-1\)來決定 (以Rule 50為例)：

1. 取出每個細胞在\(t-1\)步時的狀態\({\tt stat}={\tt array}[t-1]\)
2. 對第\(i\)個細胞而言找出它與其相鄰另兩細胞的狀態\({\tt stat}[i-1:i+2]\)
3. 匹配\({\tt stat[i-1:i+2]}\)和Rule 50並給定細胞\(i\)在第\(t\)步的狀態

不過為了執行狀態算法的第2步其實很麻煩，所以在實現上述的算法以前，我們先來看看一個在python裡比較好的方式來做這件事。

2.1 correlate函數

在numpy提供了一個叫做correlate的函數，它接受三個輸入np.correlate(arr,win, mode="same")。對於一個任意陣列 \[ {\tt arr}=[ a_1 , a_2 , a_3 , a_4 , a_5 , a_6] \tag{2} \] 和 \[ {\tt win}=[ w_1 , w_2 , w_3 ]\tag{3} \] 輸入後會得到一個與\(\tt arr\)相同維度的新陣列\({\tt arr}^\prime\) \[ {\tt arr}^\prime=[ a_1^\prime , a_2^\prime , a_3^\prime , a_4^\prime , a_5^\prime , a_6^\prime ]\tag{4} \] 而裡面的元素由關係式 \[ \begin{eqnarray*} a_{1}^{\prime} & = & w_{1}\times0+w_{2}a_{1}+w_{3}a_{2}\\ a_{2}^{\prime} & = & w_{1}a_{1}+w_{2}a_{2}+w_{3}a_{3}\\ & & \vdots\\ a_{i}^{\prime} & = & w_{1}a_{i-1}+w_{2}a_{i}+w_{3}a_{i+1} \end{eqnarray*} \] 所決定。所以透過correlate函數，我們可以得到一個元素\(a_i\)與其鄰居\(a_{i-1}\)和\(a_{i+1}\)乘上權重win的和。
　回到細胞的狀態上，如果上一步是 \[ {\tt stat}=[ 1 , 0 , 0 , 1 , 1 , 0] \] 那麼以\({\tt win}=[1 , 1 , 1 ]\)的情況下我們會得到 \[ {\tt stat}^\prime=[ 1 , 1 , 1 , 2 , 2 , 1 ] \] 如果搭配圖1的話，一個黑色代表\(1\)，兩個黑色代表\(2\)，而白色則是\(0\)，我們便可以對應到Rule 50下一步每個細胞狀態的規則。不過等等，好像哪裡怪怪的，在圖1裡，我們發現\(\tt 110\)和\(\tt 101\)都是對應到correlate的結果是\(2\)，但它們細胞在下一步的狀態完全是不一樣的，所以將權重都取成\(1\)其實不是很好的方式。
　但如果將之改為\({\tt win}=[4 , 2 , 1]\)就不一樣了，我們發現當上一步是\(\tt 000\)時，correlate對中間的元素給出\(0\)，而\(\tt 001\)時給出\(1\)，\(\tt 010\)時給出\(2\)⋯⋯到\(\tt 111\)時給出\(7\)，這些數字剛好就是圖1由右向左數的順序，所以我們便將原來要match pattern的方式改成match order就可以了，其中order是任意Rule給出的順序，這也是我們在上一個計算規則表的函數rule_tab時為什麼要加上[::-1]讓輸出的順序顛倒的緣故。我們便將狀態算法修正如下：

1. 取出每個細胞在\(t-1\)步時的狀態\({\tt stat}={\tt array}[t-1]\)
2. 使用\(\tt np.correlate\)輸出\(\tt stat\)的關聯係數 (此關聯係數便是任意Rul的順序)
3. 按關聯係數對應到\(\tt rule\_tab\)的順序 (例如對第\(i\)個細胞關聯係數是4，那麼便取出\(\tt rule\_tab\)的第4個元素的值，這個值便是第\(i\)個細胞在下一步的狀態)

使用python的腳本得到

def state(array,i,rule_table,window = [4,2,1]):
    row = array[i-1]
    corr = np.correlate(row,window,mode="same")
    # input rule_table is generated from function rule_tab
    array[i]=rule_table[corr]

這便可以計算在任意步每個細胞的狀態。

2.2 構造一個元胞自動機

有了上面兩個函數rule_table來生成任意Rule的規則表以及計算任意步的細胞狀態state後，我們便擁有了所有需要構造CA的算法了。這個CA的算法只需要兩個主要的輸入，一個是Wolfram規則Rule，介於0~255之間，另一個是要讓細胞演化幾步，CA算法可以表示如下：

1. 輸入規則\(\rm rule\)與步數\(\rm step\)
2. 依照輸入產生一個初始的零矩陣\(\tt arr\)，維度是\(({\rm step},2{\rm step}-1)\)
3. 隨機給定\({\tt arr}[0]\)每個細胞的初始狀態 (\(0\)或\(1\))，維度是\((2{\rm step}-1,)\)
   eg. \({\tt arr}[0]={\tt np.random.randint(2,size=2*{\rm step}-1)}\)
4. 計算規則\({\tt tab}={\tt rule\_tab}\)
5. 呼叫\(\tt state\)計算從第\(i=1\)步到\(i={\rm step}\)步的所有細胞狀態：
   \({\tt for}~~i~~{\tt in~~range(1,{\rm step})}:{\tt do~~state}({\tt arr},i,{\tt tab})\)
   (如果仔細看\(\tt state\)函數的話，假設做到\(i=5\)的話，\(\tt state\)會將這步得到的結果存進\({\tt arr}[5]\)裡，等於是不斷更新初始矩陣\(\tt arr\)直到第\(\rm step\)步結束為止)

上述可以用以下的腳本來實現：

def CA(rule,step):
    
    """
    CA(rule, time)
    
    Calculate the array of cellular automaton.
    
    Parameters
    ----------
    
    rule: Stephen Wolfram's rule on the cellular
    automata. It lies in between 0 and 255.
    
    time: How many time steps in the calculation.
    
    """ 
    array = np.zeros((step,2*step-1),dtype = np.uint8) 
    array[0] = np.random.randint(2, size = 2*step-1)
    tab = rule_tab(rule)
        
    for i in range(1,step):
        state(array,i,tab)

    return array

最後我們可以用imshow來繪製CA的結果，例如

# Calculating Rule 129 with 201 steps  
ca129 = CA(129,201)
    
# Plot
plt.imshow(ca129,cmap = "binary",interpolation = "None")

結果如圖3所示，在圖3裡我們特別讓初始條件是只有中間那個細胞的值是\(1\)其它都是\(0\)開始。

圖3. Rule 119產生201步的CA，初始條件為中心位置的細胞狀態為\(1\)其餘皆\(0\)

　在github的notebook裡，我在函數CA還放入自動繪圖和初始條件是隨機或是永遠保持中間的細胞狀態為\(1\)其餘是\(0\)的選擇，有興趣可以自己看。題外話是在知名的商用數學套件Mathematica®裡有一條函式CellularAutomaton可以生成任意規則的CA以及RulePlot來視覺化不同的Rule，如圖1和圖2等，有興趣可以詳見該套件的參考資料。

三、回顧元胞自動機的分類

這個小節我們很簡單的回顧一下Stephen Wolfram對CA的4個分類，我們先看下放用4個不同Rule所產生201步的CA。

圖4. 不同Rule產生的CA

這4張圖的初始條件都是隨機給定的，亦即細胞在\({\tt arr}[0]\)可以任意處在\(0\)或\(1\)的狀態上。Rule 129和110沒有意外是屬於第3類的狀態，原來的隨機性消失了，取而代之是其他的隨機噪聲，沒有任何明顯的結構長存，或至少在201步的步伐是看不出來週期性的結構。而Rule 210是屬於第2類，圖像右半部的黑白斜線是一開始初始條件給定的隨機性，在往後的演化中被保留下來的，但是發生了位移，而那些在右邊界消失的初始隨機性的圖騰則逐漸被CA本身的三角形結構所取代。那你覺得Rule 20和圖3的Rule 119是屬於哪一類呢？

四、小結

CA算是已被提出近一個世紀的概念了，雖然目前尚未在物理世界中有直接顯著的對應例子，但即使是這種純粹的柏拉圖世界的圖像，我們也能從很簡單的幾個規則出發，逐漸構造出千變萬化又令人著迷的圖像！


Source on github: I&A5_cellular_automaton.ipynb