Agenda Cluster Monte Carlo

Main.hs

module Main where

import qualified Data.Vector.Unboxed as V
import qualified Data.Vector.Unboxed.Mutable as M

import Control.Monad.Random

import System.Environment

import Text.Printf

shuffle vector = V.thaw vector >>= shuffle' (V.length vector)
  where
    shuffle' i v | i <= 1 = V.unsafeFreeze v
                 | otherwise = do
                     let i' = i - 1
                     index <- getRandomR (0,i')
                     M.swap v i' index
                     shuffle' i' v
        
makeDeck cards agendas = V.fromList [x <= agendas | x <- [1..cards]]

checkClump clumpSize agendas deck = go 0 0
  where
    go i n | i >= V.length deck = False
           | otherwise = let n' = if deck V.! i then n+1 else n
                             n'' = if i >= clumpSize && deck V.! (i-clumpSize) then n'-1 else n'
                         in if n'' >= agendas
                            then True
                            else go (i+1) n''

runTrials successes trialsLeft _ _ _ | trialsLeft < 1 = return successes
runTrials successes trialsLeft clumpSize agendas deck = do
  deck' <- shuffle deck
  let succ' = if checkClump clumpSize agendas deck' then successes + 1 else successes
  runTrials succ' (trialsLeft-1) clumpSize agendas deck'
  
main :: IO ()
main = do
  args <- getArgs
  let deckSize = read (args !! 0) :: Int
      agendaCount = read (args !! 1) :: Int
      clumpSize = read (args !! 2) :: Int
      clumpAgendas = read (args !! 3) :: Int
      trials = read (args !! 4) :: Int
      
  succ <- runTrials 0 trials clumpSize clumpAgendas (makeDeck deckSize agendaCount) :: IO Int

  printf "Deck: %d/%d\nOver %d trials we had a %d/%d agenda clump %d times (%.1f%%)\n"
    deckSize agendaCount trials clumpAgendas clumpSize succ (100 * fromIntegral succ / fromIntegral trials :: Double)