This week only. Pastebin PRO Accounts Christmas Special! Don't miss out!Want more features on Pastebin? Sign Up, it's FREE!
Guest

Ising Metropolis in Haskell

By: a guest on Mar 7th, 2010  |  syntax: Haskell  |  size: 3.34 KB  |  views: 810  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. module Main where
  2.  
  3. import System.Random.Mersenne as M
  4. import Data.Vector.Unboxed.Mutable as U
  5. import Control.Monad.State as S
  6. import Control.Monad.Primitive
  7. import Control.Applicative
  8. import Control.Concurrent
  9. import Graphics.UI.SDL as SDL
  10.  
  11. newtype Grid = Grid { unG :: (Int, MVector (PrimState IO) Int) }
  12. newtype Ising = Ising { unI :: (M.MTGen, Grid) }
  13.  
  14. type Pos = (Int, Int)
  15.  
  16. (|>) :: Grid -> Pos -> IO Int
  17. (Grid (n,g)) |> (x,y) | x == n     = U.unsafeRead g (y*n)
  18.                       | x == -1    = U.unsafeRead g (y*n + n - 1)
  19.                       | y == n     = U.unsafeRead g x
  20.                       | y == -1    = U.unsafeRead g ((n - 1)*n + x)
  21.                       | otherwise  = U.unsafeRead g (y*n + x)
  22.  
  23. (<==) :: Pos -> Int -> (Pos, Int)
  24. pos <== val = (pos,val)
  25.  
  26. (<|) :: Grid -> (Pos, Int) -> IO ()
  27. Grid (n,g) <| ((x,y), val) = U.unsafeWrite g (y*n + x) val
  28.  
  29. infixr 1 <==
  30. infixr 0 <|
  31.  
  32. -- Just a shortcut
  33.  
  34. l = lift
  35.  
  36. -- Initialize a ising system
  37.  
  38. initIsing :: M.MTGen -> Int -> IO Ising
  39. initIsing rg n = do
  40.             g <- U.new $ n*n
  41.             grid <- return $ Grid (n, g)
  42.             let spin x = x + (x - 1)
  43.             sequence_ [ do rnumber <- random rg
  44.                            grid <| (i,j) <== spin (rnumber `mod` 2) | i <- [0..n - 1], j <- [0..n - 1] ]
  45.             return $ Ising (rg, grid)
  46.  
  47. -- Update the grid of Ising Model 'ns' steps in temperatur 't'
  48.  
  49. updateIsing :: Int -> Double -> StateT Ising IO ()
  50. updateIsing ns t = up 0 ns
  51.     where
  52.           up i ns | i == ns = return ()
  53.                   | otherwise = do
  54.                         Ising (rg, grid@(Grid (n,_))) <- S.get
  55.                         x <- l $ (`mod` n) <$> random rg
  56.                         y <- l $ (`mod` n) <$> random rg
  57.                         (e,s) <- l $ calc_energy grid (x,y)
  58.                         if e < 0
  59.                             then do l $ grid <| (x,y) <== negate s
  60.                             else do p <- l $ (random rg :: IO Double)
  61.                                     when (p < exp (-2*(fromIntegral e)/t)) $ do
  62.                                         l $ grid <| (x,y) <== negate s
  63.                         up (i + 1) ns
  64.  
  65.           calc_energy g (x,y) = do
  66.               s <- g |> (x,y)
  67.               e <- (s*).(sum) <$> (mapM (g |>) [(x, y - 1), (x, y + 1), (x - 1, y), (x + 1, y)])
  68.               return (e, s)
  69.  
  70. -- Draw ising
  71.  
  72. drawIsing screen (Ising (_, grid@(Grid (n,_)))) = do
  73.     sequence [ drawSquare screen i j | i <- [0..n - 1], j <- [0..n - 1] ]
  74.     SDL.flip screen
  75.         where drawSquare s i j = do
  76.                 s <- grid |> (i,j)
  77.                 let color 1 = 0xFF0000
  78.                     color _ = 0x00FF00
  79.                     spSize = screen_size `div` n
  80.                     rect i j = Just $ Rect (i*spSize) (j*spSize) spSize spSize
  81.                 SDL.fillRect screen (rect i j) (SDL.Pixel $ color s)
  82.                 return ()
  83.  
  84. -- Update ising and draw
  85.  
  86. sdlUpIsing s i t c = do
  87.     evalStateT (updateIsing 1000 t) i
  88.     drawIsing s i
  89.     threadDelay 1
  90.     sdlUpIsing s i (t - 0.001) (c + 1)
  91.  
  92. screen_size = 400
  93.  
  94. main = do
  95.     SDL.init [InitEverything]
  96.     setVideoMode screen_size screen_size 32 []
  97.     screen <- SDL.getVideoSurface
  98.     rg <- M.newMTGen Nothing
  99.     i <- initIsing rg 200
  100.     forkIO . forever $ waitEvent >>= \e -> when (e == Quit) quit
  101.     sdlUpIsing screen i 3 0
clone this paste RAW Paste Data