IA6 - 小專題: 進擊的黴菌

初始3個菌點經過80單位時間後形成的菌落

利用電腦模擬真實世界的情況，不僅僅是電腦科學家，還包含了物理學家、生物學家、化學家和多種領域專家有志一同的想法。自然界千變萬化，圖騰多樣，但或許它們都可以從幾個簡單的規則當中得到答案，這便是所謂的湧現現象 (emergence phenomenon)。我們來看個簡單的例子，黴菌菌落的增值過程。

一、菌落的消長

用一個\(N\times N\)的區域代表一整片的培養區域如圖1，這個區域裡有比例\(q\)的部份是培養基質可以生長出黴菌 (淺綠)，剩下的區域則不能長出黴菌 (白色)，假設長出和不能長出的區域是隨機分佈的，那麼我隨意在這片區域灑上幾個菌點 (墨綠色)，是問它們發展成菌落的過程為何？

圖1. 培養區域，淺綠是含有生長基質的區域，白色不能生長，墨綠色表示初始的5個隨機菌點

　這個問題和Prof. Downey在他的著作Think Complexity的第七章在描述滲透的情況基本上一模一樣，除了關鍵的黴菌生長需要引入額外的判定條件外，所以我會依照Prof. Downey的方式來做，最後在引入這個判定生長的條件。欲執行這個雜記的notebook，我需要下載Prof. Downey所寫的一個python腳本Cell2D.py，可以從雜記的notebook裡直接點連結。這個腳本主要是需要其中的Cell2DViewer類來將黴菌菌落的消長輸出成動畫，如果只需要用imshow看特定時刻的結果，引入Cell2D.py不是必須的，但如果想執行本篇雜記的notebook第一個擴散的例子，則需要有其中的Cell2D類。

1.1 預備知識

如果說以圖1其中一個淺綠的基質為中心 (Eq. (1)的\(?\)處)，那麼它週邊\(3\times 3\)的區域內4個東南西北 (不含對角線方向：東南、東北……等等) 的方向可能含有黴菌或不含有黴菌，如下 \[ \begin{array}{ccc} \times & {\rm 黴} & \times\\ {\rm 無} & ? & {\rm 無}\\ \times & {\rm 無} & \times \end{array}\tag{1} \] 所以我們想知道如果它周圍有一個以上的菌點，那麼它會不會被侵入也長出黴菌？
　若有時間閱讀Think Complexity的第七章的話，這原則上是一個滲透的問題，只要乾燥處的週邊有一個溼掉的點，那麼乾燥處經過一個單位時間後就一定會被滲透弄溼，機率是\(100%\)。不過黴菌的生長則不然，就算這個區域含有黴菌生長所需要的基質，鄰近的區域剛好也有一個黴菌，也不代表這個區域在下一個單位時間內就一定會長出黴菌，而應個要用機率的方式來決定。就像是班上有人感冒，我坐在感冒同學旁邊的不一定馬上就會被感染一樣，會不會被感染是機率，當然旁邊有一個或愈多個感冒的同學，我得病的機率就愈大。
　回到黴菌的問題上，假設在下一個單位時間後中心區域會長出黴菌我們賦與它一個值\(1\)而不長出黴菌則是\(0\)，我們可以採用隨機抽樣的方式來決定，抽到\(1\)則長出抽到\(0\)則沒事。那麼抽到決定長出或不長出黴菌的機率就用該中心區域的週邊有多少黴菌來給定。因為東南西北共4個方向，如果有一個菌點存在則抽到\(1\)的機率是\(1/4\)，兩個是\(1/2\)以此類推，我們變可以用

np.random.choice([1,0], p = [i/4, 1-i/4])

來決定該處是否在下一刻長出黴菌，其中\(i\)是鄰近區域有幾個菌點，至多4個。
　而這個判別便是與Think Complexity第七章的滲透所不一樣的地方，再來滲透問題的起始條件是固定的，原則上是從最上排開始往下滲 (大家都無法抗拒地心引力)，但黴菌菌落的生長是從隨機選取的區域開始！我們稍候會再回到如何寫一個判別被侵入的函數的問題上。

二、建立黴菌的環境

一開始我們先建立一個黴菌的類 (class) 叫作mold，而初始化是形成一個\(N\times N\)的陣列，其中陣列的元素可以在\(0\)和\(1\)兩者隨機取值，是\(1\)代表可長出黴菌的基質，\(0\)代表不能長出東西的地方。

class mold():
    
    window = np.array([[0, 1, 0],
                       [1, 0, 1],
                       [0, 1, 0]])
        
    def __init__(self,n,q,m):        
        self.q = q
        self.array = np.random.choice([1,0],(n,n),p=[q,1-q])

        # Randomly put mold seeds.
        # If the location has element 5, then there exists a mold. 
        for i in range(m):
            self.array[np.random.randint(n),np.random.randint(n)]=5

在__init__裡初始化的陣列中可讓黴菌生長的基質佔整個區域的比例q，而剩下長不出黴菌的則佔了1-q，最後一步則是隨機灑m個菌點。所以對mold初始化便是建立黴菌生長的區域，以及菌點的數量。出來的結果應該跟圖1相仿，它是一個\(20\times 20\)的區域，其中基質佔q=0.8而初始菌點m=5。另外因為\(1\)和\(0\)已經用來描述可不可以長出黴菌的代表，我們讓長出黴菌的地方用\(5\)來代稱。
　當環境建立好後，我們變可以計算下一刻後從這個菌點發展出去的菌落。按照上一節的作法，我們來以個陣列元素為中心，像Eq. (1)一樣，我們可以去計算這個中心的週邊4個方向的情況，一定是0、1和5三種擇一。這可以透過scipy裡的correlate2d來完成，基本上它就和np.correlate做一模一樣的事，但現在變成二維的形式，見圖2。

圖2. correlate2d示意圖

　給定window的大小，例如圖2取\(3\times 3\)，所以以中間為核，我要看它鄰近的8個元素依權重給定的關聯值。圖二右當kernel掃到對應的中間是元素2，而這個2的關聯值便由它周圍的六個元素依對應的權重值\(w_i\)給定，故對2而言有 \[ {\rm correlation}=5w_1+w_2+7w_3+3w_4+5w_2+9w_6+4w_7+7w_8. \] 　不過對我們要計算鄰近的菌點數不需要對角方向，所以在mold類裡window只有東南西北4個方向的權重\(1\)，其他都是\(0\)。當週邊有一個菌點以上的話，關連函數會給出該中心對應的關連值大於等於\(5\)，而4個方向都存在菌點有最大是\(20\)，沒有菌點存在的話則小於\(5\)。

2.1 構造侵入函數

當決定周邊有有菌點的時候，我們就計算中間區域被該菌點侵入的機率，回到上一節繼續我們的侵入函數。我先用correlate2d得出當下整個區域每個陣列元素對應的關聯值，這個關聯值不僅僅告訴我們他旁邊有幾個菌點，這個值還告訴我們他被侵入的機率。由於我知道周邊都被黴菌填滿的話有最大值\(20\)，所以若關聯值是\(20\)，則我讓該區域被侵入的機率是百分百，若小於\(20\)比如說是11 (代表有兩個菌點、一個可長出黴菌的基質和一個長不出黴菌的點：\(5\times 2+1\times 1+1\times 0 =11\))，則被侵入的機率是\(11/20\)。所以我可以自訂侵入函數penetration，其輸入是整個關聯陣列，而輸出的侵入陣列是與關聯陣列同維度，每個元素都是\(1\)和\(0\)，分別代表該處在下一步會被侵入或不被侵入：

# Define penetration function
def penetration(corr):    
    row,col = corr.shape    
    pen = np.zeros((row,col))     # Initializing an empty penetration matrix
    
    for r in range(row):
        for c in range(col):
            pen[r,c] = np.random.choice([1,0],p=[corr[r,c]/20,1-corr[r,c]/20])
 
    return pen

在penetration函數內的for迴圈檢查的是關聯陣列每個位置的關聯值corr[r,c]，再用這個關聯值當作被侵入機率corr[r,c]/20。

2.2 完成黴菌的演化

接著我們要繼續完成mold類，由於我們只定義了初始化後環境的狀態 (哪裡可以/不可以長黴以及隨機撒的菌點)，我們接著要定義每一步後黴菌在這個區域的狀態的步驟函數step。
　首先先找出關聯陣列，再利用關聯陣列取得侵入陣列，最後和初始陣列做比較 (這裡的初始不見得是mold類最初初始化後得出那個只有幾個菌點的陣列，也可以是經過多步演化後得到的狀態陣列)，當這三個陣列的同一位置滿足下面三個條件時：

1. 該處目前沒有黴菌，初始陣列上的元素是\(1\) (如果是\(0\)的話長不出黴菌)
2. 該處的關聯值必須大於等於\(5\) (週邊至少一處有菌點)
3. 侵入陣列在該處的值為\(1\) (代表它會被鄰近的菌點侵入)

當這3個條件都滿足後，我們便將初始陣列該處的值更新為\(5\)，代表他長出黴菌了，否則就維持原來的值不動。所以步驟函數step可以寫成

# Calculating the evolving of mold at given status
    def step(self):           
        a = self.array
        corr = correlate2d(a,self.window,mode='same')      # finding correlation
        pen = penetration(corr)                            # finding penetration

        # checking if the 3 conditions are satisfied altogether
        # if so, updates the value by 5 at that location
        # if not, it remains untouched
        self.array[(a==1) & (corr>=5) & (pen==1)] = 5

所以我們便完成mold類的撰寫了，現在我們可以用它來計算給定菌落的狀態後下一步演化的結果。

三、黴菌，始動！

該建的都建完了，模擬就簡單啦！先初始化一個長黴菌的區域mm，有\(65\%\)的區域覆蓋可供黴菌生長的基質，一開始的菌點為\(5\)個：

mm = mold(n=50,q=0.65,m=5)
plt.imshow(mm.status(),cmap="Greens",vmin=0,vmax=5)

得到了圖2。上面我是在mold類裡多定義一個status函數，它本身沒什麼用處，就僅僅是回傳當前這個區域陣列的狀態而已，由於當前的狀態是放在該類底下的變數self.array裡，所以也可以直接輸入mm.array存取，跟mm.status()得結果完全一致，只是用mm.status()看起來比較假掰而已😂

圖3. 一開始mm的狀態

　當我們有一開始的mm後，我們便可以利用mm.step()來更新演化到下一步時mm的狀態，假設我要看演化到第200步的狀態，我可以用

for i in range(201):
    mm.step()

再用imshow檢視mm.status()即可。圖4展示了mold(n=200,q=0.65,m=10)演化到不同步數時的狀態。

圖4. 不同時間時黴菌菌落增殖的情況

　另外我們也可以用Prof. Downey的Cell2DViewer製作動畫，結果如圖5，怎麼做請直接看notebook。

圖5. mold(n=200,q=0.65,m=5)隨時間演化的動畫

四、小結

這裡我基於Think Complexity第七章滲透的內容製作出黴菌增長的模擬，差異在於Prof. Downey在討論滲透時不需要去知道有沒有被侵入的機率，亦即不需要引入penetration函數。這是可以理解的，因為只要鄰近被浸濕的話，那麼不可避免的該處必然會濕。但黴菌的增長則不然，即便周遭發霉了，也不代表中心處就會發霉，我僅能用機率的方式說它發霉的機率比較高，特別是周遭發霉的數量愈多的話，而這也是我們這篇雜記和Think Complexity第七章滲透的主題不一樣的地方。如果把penetration函數拿掉的話，會發現黴菌菌落增值的速度快很多，why?
　但我必須強調即使用penetration的方式也仍然只是一種近似，真實的菌落不僅僅是增而已，還有可能會減，甚至當菌落數變大之後，他能攜帶孢子散佈出去，使得演化到後期會有一些地方本來沒有發黴卻突然冒出菌點，這都是這裡所沒有考慮的，當然更有多種菌落彼此競爭的情況，這都會使模擬變得複雜，但我相信也會變得更有趣！


Source on github: I&A6_mold.ipynb