用 Mathematica 搜索生命游戲中的靜物（二）

之前寫的那篇《用 Mathematica 搜索生命游戲中的靜物》谜嫉，由于自己寫的部分多了一點（雖然關鍵的一步還是用 Mathematica 自帶的 FindPath 函數(shù)），特別慢谓形，還特別耗內存泡躯。果然對我這種完全不懂算法的人，就應該把所有的事情都交給 Mathematica 才對阱佛。

生命游戲里的每個細胞的狀態(tài)可以看成一個布爾值台丛，一個（大小有限的）圖樣則可以看成一個由布爾值組成的二維數(shù)組耍缴。于是，要尋找滿足某些條件圖樣挽霉，就相當于要求這些布爾值滿足某個方程防嗡。于是，這是一個布爾可滿足性問題（SAT）侠坎。Mathematica 里有個叫 SatisfiabilityInstances 的函數(shù)就是干這個的蚁趁。

我們先把要滿足的條件用 Mathematica 代碼寫出來。

假設我們要找的圖樣在一個 x 乘 y 的長方形里邊实胸，它對應的數(shù)組可以用 Array[b, {x, y}] 來表示他嫡，這里每個 b[i,j] 表示一個細胞。

靜物（Still life）指的是穩(wěn)定的圖樣庐完，它要滿足下面兩個條件：

每個活細胞周圍的活細胞個數(shù)必須是2或者3钢属，
每個死細胞周圍的活細胞個數(shù)不能是3。

我們先寫一個函數(shù)來數(shù)一個細胞周圍的活細胞個數(shù)：

NeighborCount[k_, {i_, j_}] :=
  BooleanCountingFunction[{k},
   Delete[Catenate[Array[b, {3, 3}, {i, j} - 1]], 5]];

然后每個細胞要滿足的條件可以寫成這樣：

StillLifeCondition[i_, j_] :=
  (b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
   (! b[i, j] && ! NeighborCount[{3}, {i, j}]);

此外门躯，這個圖樣是有限的淆党，邊界之外的細胞總是死的。于是讶凉，整個圖樣要滿足的條件可以寫成：

Array[StillLifeCondition[##] /.
   b[i_, j_] /; i < 1 || i > x || j < 1 || j > y :> False &,
 {x, y} + 2, 0, And]

然后我們可以用 SatisfiabilityInstances 函數(shù)：

SearchStillLife[x_, y_] :=
 ArrayReshape[
  SatisfiabilityInstances[
    Array[StillLifeCondition[##] /.
       b[i_, j_] /; i < 1 || i > x || j < 1 || j > y :> False &,
     {x, y} + 2, 0, And],
    Catenate[Array[b, {x, y}]]][[1]], {x, y}];

然后我們就可以用這個 SearchStillLife 函數(shù)來找靜物染乌。不過，它出奇地慢懂讯，找個16乘16大小的靜物也要花上44秒：

ArrayPlot[Boole@SearchStillLife[16, 16], Mesh -> All] // AbsoluteTiming

能不能快一些荷憋？

我不清楚這個 SatisfiabilityInstances 函數(shù)用的是什么算法。不過看了一下維基百科褐望，SAT 問題最常用的好像是一種叫 DPLL 的算法勒庄。我完全不懂算法串前，就不管它的細節(jié)了」Γ總之，它要求輸入的是一個合取范式（CNF）减宣。

于是盐须，我們可以試試把條件轉換成 CNF，也就是說漆腌，把前面的 StillLifeCondition 函數(shù)改寫成：

StillLifeCondition[i_, j_] :=
 BooleanConvert[
  (b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
   (! b[i, j] && ! NeighborCount[{3}, {i, j}]), "CNF"];

后面的 SearchStillLife 函數(shù)定義不變≡舻耍現(xiàn)在再找一遍16乘16的靜物，果然快了很多闷尿，只花了2.2秒塑径。不過找出來的靜物……

居然是空的……

看來 SatisfiabilityInstances 函數(shù)在處理 CNF 時用的是和別的形式不同的算法。滿足條件的靜物有很多填具，但 SatisfiabilityInstances 只會輸出其中的第一個统舀。而在輸入 CNF 的時候，不巧空靜物就是這“第一個”劳景。

只能找到空靜物的代碼并沒有什么用誉简。我們可以試試給原來的代碼引入一點隨機的因素，讓這個空靜物不再是“第一個”盟广。比如說闷串，把原來數(shù)組異或上一個隨機的數(shù)組，把前面的 SearchStillLife 函數(shù)改寫成：

SearchStillLife[x_, y_] := 
 Block[{r = RandomChoice[{True, False}, {x, y}]},
  MapThread[Xor, {r,
    ArrayReshape[
     SatisfiabilityInstances[
       Array[StillLifeCondition[##] /. 
           b[i_, j_] /; i < 1 || i > x || j < 1 || j > y :> False /.
                   b[i_, j_] :> Xor[b[i, j], r[[i, j]]] &,
        {x, y} + 2, 0, And],
       Catenate[Array[b, {x, y}]]][[1]], {x, y}]}, 2]]

現(xiàn)在完整的代碼變成了：

NeighborCount[k_, {i_, j_}] :=
  BooleanCountingFunction[{k},
   Delete[Catenate[Array[b, {3, 3}, {i, j} - 1]], 5]];

StillLifeCondition[i_, j_] :=
 BooleanConvert[
  (b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
   (! b[i, j] && ! NeighborCount[{3}, {i, j}]), "CNF"];

SearchStillLife[x_, y_] := 
 Block[{r = RandomChoice[{True, False}, {x, y}]},
  MapThread[Xor, {r,
    ArrayReshape[
     SatisfiabilityInstances[
       Array[StillLifeCondition[##] /. 
           b[i_, j_] /; i < 1 || i > x || j < 1 || j > y :> False /.
                   b[i_, j_] :> Xor[b[i, j], r[[i, j]]] &,
        {w, h} + 2, 0, And],
       Catenate[Array[b, {w, y}]]][[1]], {x, y}]}, 2]]

這次能找到好看的靜物了筋量，花的時間時間在2.7秒左右：

SeedRandom[233];
ArrayPlot[Boole@SearchStillLife[16, 16], Mesh -> All] // AbsoluteTiming

再試試大一點的靜物烹吵，比如說64乘64的〗拔洌花了大概46秒肋拔。

找完了靜物，再來找振蕩子（Oscillator）呀酸。振蕩子是隨時間周期變化的圖樣只损。于是，我們可以給那個布爾值的數(shù)組增加一個時間的維度七咧，寫成 Array[b, {p, w, h}]跃惫，這里 p 是它的周期。它們滿足的條件也要作出相應的修改：

NeighborCount[k_, {t_, i_, j_}] := 
  BooleanCountingFunction[{k}, 
   Delete[Flatten[Array[b, {1, 3, 3}, {t, i - 1, j - 1}]], 5]];

OscillatorCondition[t_, i_, j_] := 
  BooleanConvert[
   (b[t, i, j] && (b[t + 1, i, j] ? 
        NeighborCount[{2, 3}, {t, i, j}])) ||
    (! b[t, i, j] && (b[t + 1, i, j] ? 
        NeighborCount[{3}, {t, i, j}])), "CNF"];

SearchOscillator[p_, w_, h_] :=
 Block[{r = RandomChoice[{True, False}, {p, w, h}]},
  MapThread[
   Xor, {r, 
    ArrayReshape[
     SatisfiabilityInstances[
       Array[OscillatorCondition[##] /.
           {b[t_, i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False, 
            b[0, i_, j_] :> b[p, i, j]} /.
          b[t_, i_, j_] :> Xor[b[t, i, j], r[[t, i, j]]] &,
        {p, w + 2, h + 2}, 0, And],
       Flatten[Array[b, {p, w, h}]]][[1]], {p, w, h}]}, 3]]

（代碼里的這個 ? 符號表示的是等價艾栋，Mathematica 里顯示為 ?爆存。）

可能因為靜物比較稀少，速度比找靜物時慢了很多蝗砾，找一個周期2的16乘16的靜物花了11秒：

SeedRandom[233];
ArrayPlot[#, Mesh -> All] & /@ 
  Boole@SearchOscillator[2, 16, 16] // AbsoluteTiming

周期3的就更慢了先较，找下面這個振蕩子花了14分鐘携冤。

更高的周期我就不敢試了。

應該會有更快的辦法闲勺。有好的建議歡迎評論曾棕。