読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

モンテカルロ法

コンピュータ囲碁モンテカルロツリー探索に興味津々でしたが、コーディングがわからずほったらかしでした。そもそもモンテカルロ法さえ組んだことがない。
西尾さんという方のブログの記事「Rediscover the Monte Carlo」を読ませていただいて、これならHaskellに焼き直しできるかもと、やってみました。

import MyBase ((|$))
import System.Random (getStdGen, randomRs)
import Data.Array
import Data.List (insertBy)
import Data.Ord (comparing)

trialsPerOutput = 100000
timesOfOutputs = 20

splitPer n xs = let (xs1, xs2) = splitAt n xs in xs1 : splitPer n xs2

monte rs = take timesOfOutputs $ zipWith (/) (scanl1 (+) averages) [1..]
  where
    average ns = sum ns / fromIntegral (length ns)
    averages = rs |$ map (^2) |$ splitPer 2 |$ map sum |$ 
               map (\x -> if x < 1.0 then 4.0 else 0.0) |$ 
               splitPer trialsPerOutput |$ map average

score a b =  let 
  a1 = fromIntegral a + 1
  b1 = fromIntegral b + 1
  n = a1 + b1 + 2
  in
   (2 * a1 * b1 + n) / n / (n + 1) / (n + 2)

width = 100

oneTurn (trials, hits, ((_, (ix, iy)) : ps)) r1 r2 = let
  nextTrials = accum (+) trials [((ix, iy), 1)]
  x = r1 / fromIntegral width + 1 / fromIntegral width * fromIntegral ix
  y = r2 / fromIntegral width + 1 / fromIntegral width * fromIntegral iy
  nextHits = if x^2 + y^2 < 1 then accum (+) hits [((ix, iy), 1)] else hits
  s = score (nextHits ! (ix, iy)) (nextTrials ! (ix, iy) - nextHits ! (ix, iy))
  priorities = insertBy (flip (comparing fst)) (s, (ix, iy)) ps
  in
   (nextTrials, nextHits, priorities)
  
turns init rs = aux init rs
  where
    aux state (r1 : r2 : rs) = let 
      next = oneTurn state r1 r2 
      in
       next : aux next rs
  
area trials hits = 4 * (sum $ zipWith (/) (map fromIntegral (elems hits)) (map fromIntegral (elems trials))) / (fromIntegral width ^ 2)
  
monte2 rs = let
  zeros = array ((0, 0), (width - 1, width - 1)) [((x, y), 0) | x <- [0..width - 1], y <-[0..width - 1]]
  initPriorities = [(score 0 0, (x, y)) | x <- [0..width - 1], y <- [0..width - 1]]
  in
   turns (zeros, zeros, initPriorities) rs |$ thinOut trialsPerOutput |$ 
   map (\(t, h, _) -> area t h) |$ take timesOfOutputs
  where
    thinOut n xs = let (pres, posts) = splitAt n xs in head posts : thinOut n posts

main = do
  g1 <- getStdGen
  g2 <- getStdGen
  let rs1 = (randomRs (0.0, 1.0) g1) :: [Float]
  let rs2 = (randomRs (0.0, 1.0) g2) :: [Float]
  print $ monte rs1 
  print $ monte2 rs2  

monte2は大きなArrayの更新を繰り返すので、コピー繰り返してだめだろうなー、でももしコンパイラが賢かったら、内部で破壊的に更新してくれてすごい速いかも、とちょっと期待しながらやってみたら、やっぱりダメでした…2GHz Core 2 Duoで40分かかります。
STモナドですか。副作用を閉じ込めるためのモナドなら少しは納得がいきますが、参照透過なコードで書けるところを、実行速度を上げるためにモナド使うのって、それなら最初から代入ありの言語使うよなー。
Haskellにこだわるか、Haskellから離れるか微妙な心境です。