起點

以下描述一個需要 heavy optimisation 的 Haskell 模擬程式。這是 Weijin 最近加入中研院物理所李副所長旗下接到的第一個習題，用電腦模擬某種 Ising model 並作圖。對我而言這是練習 practical Haskell 的大好機會 XD。

一開始定義「晶格」為 n * n 的矩陣，每一格的值可能是 +1 或 -1。這裡直接用 array of Bools 表示。

type Lattice = IOUArray (Int, Int) Bool

可以把 lattice 印成 0 和 1 的矩陣。

showLattice lat n = getElems lat >>= putStr . concat .
  zipWith (flip (++)) (cycle $ replicate (n-1) "" ++ ["\n"]) .
  map (show . fromEnum)

給定一個 lattice，可以計算某一格的能量：把這一格的值和它上下左右四個鄰居分別相乘後加起來取負值。Lattice 的左右、上下邊界是黏在一起的，例如左邊界上的格子的左邊鄰居就是同一列右邊界的格子。物理上叫它作 "periodic boundary condition"，拓樸上這個 lattice 就變成一個甜甜圈，日常生活的例子就是不撞牆的貪食蛇 XD。下面這個函式取得某個座標和它的鄰居的值。

neighbors :: Lattice -> Int -> (Int, Int) ->
             IO (Bool, Bool, Bool, Bool, Bool)
neighbors lat n (i, j) =
  do c <- readArray lat (i, j)
     let dec x = let xdec = x - 1 in if xdec <  0 then n - 1 else xdec
         inc x = let xinc = x + 1 in if xinc >= n then 0     else xinc
     u <- readArray lat (dec i, j)
     d <- readArray lat (inc i, j)
     l <- readArray lat (i, dec j)
     r <- readArray lat (i, inc j)
     return (c, u, d, l, r)

然後是計算能量的函式。

energy :: Lattice -> Int -> (Int, Int) -> IO Int
energy lat n (i, j) =
  do (c, u, d, l, r) <- neighbors lat n (i, j)
     let f b1 b2 = if b1 == b2 then 1 else -1
     return $! -(f c u + f c d + f c l + f c r)

可以印出整個 lattice 的平均能量。

showAverage lat n =
  do r <- getBounds lat
     xs <- sequence (map (energy lat n) (range r))
     print (fromIntegral (sum xs) / fromIntegral (n^2))

這個模擬程式主要要做的事情是這樣：隨機選一格，算出它目前的能量和翻轉它（1 和 -1 互換）之後的能量。如果能量差小於（等於）0 就立刻翻轉；如果能量差大於 0 就以某個機率翻轉它。這個機率由某個溫度參數 t 和能量差 e 決定。

flipProb :: Double -> Double -> Double
flipProb t e = exp (-e / t)

下面的函式實作上述行為。

touch :: Lattice -> Int -> Double -> IO ()
touch lat n t =
  do (i, j) <- newStdGen >>= return . choose n
     (nv, e) <- energyDiff lat n (i, j)
     let flip = writeArray lat (i, j) nv
     if e <= 0
        then flip
        else do f <- newStdGen >>= return . fst . random
                if f < flipProb t (fromIntegral e)
                   then flip else return ()

（其實我不知道可不可以這樣仰賴 newStdGen。）其中 energyDiff 是 energy 的簡單改版。

energyDiff :: Lattice -> Int -> (Int, Int) -> IO (Bool, Int)
energyDiff lat n (i, j) =
  do (c, u, d, l, r) <- neighbors lat n (i, j)
     let f b1 b2 = if b1 == b2 then 1 else -1
         g c = -(f c u + f c d + f c l + f c r)
     return $! let nc = not c in (nc, g nc - g c)

choose 則是一個 pure function。

choose :: RandomGen g => Int -> g -> (Int, Int)
choose n g = let (x, _) = randomR (0, n^2 - 1) g
             in  (x `div` n, x `mod` n)

一個計算單位是 n^2 次 touch，也就是平均每格碰一次。聽說要從 n = 4 算到 n = 64，每個 n 各做兩億個單位，所以計算量驚人… 總之 main 先讀入 n、溫度、和要模擬的單位數，亂數初始化 lattice，每做完一個單位就輸出目前的平均能量，最後輸出 lattice 的狀態。

main = do n <- readLn
          t <- readLn
          cnt <- readLn
          lat <-
            newStdGen >>= newListArray ((0, 0), (n-1, n-1)) . randoms
          sequence_ $ flip map [1..cnt]
            (\i -> do replicateM_ (n^2) (touch lat n t)
                      putStr (show i ++ "\t")
                      showAverage lat n)
          showLattice lat n

目前大概慢 Weijin 的 C++ 版本至少幾百倍甚至千倍，所以這也是探索 Haskell 程式最佳化的大好機會 XD。

--
另外想寫一篇 regex 與基本 computability theory 的通識文章 XD。