$\newcommand{\defeq}{\mathrel{\mathop:}=}$

2008/08/02

起點

type Lattice = IOUArray (Int, Int) Bool


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)


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


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。）其中 energyDiffenergy 的簡單改版。
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


--

Labels: