--
-- Metropolis-Hastings Algorithm on the 1-D Ising model
-- Code by NOTtheMessiah
--
-- Parameters set in lines 28 and 54-56
--

import System.Random

--actually the negative hamiltonian for convenience
ham :: Double -> Double -> [Int] -> Double
ham j h [] = 0.0
ham j h [x] = h * fromIntegral x
ham j h (x:xs) = h * fromIntegral x + j* fromIntegral x * fromIntegral (head xs) + ham j h xs

swap :: Int -> [Int] -> [Int]
swap n a = x ++(-c : y)
    where
        (x,c:y) = splitAt n a

iff :: Bool -> a -> a -> a
iff True  x y = x
iff False x y = y

step :: [Int] -> StdGen -> ([Int],StdGen)
step a g= (better a b, g'')
    where
        (j,h,beta) = (0.025, 1.00, 1.00)
        hi = ham j h a
        hf = ham j h b
        b = swap n a
        (n,g') = randomR (0, (length a)-1) g
        p = min 1.0 $exp $ beta * (hf - hi)
        (p',g'') = randomR (0.0,1.0) g'
        better a b = iff (p' < p) a b

f = uncurry step

--initializes system
initSpin :: Int -> StdGen -> ([Int],StdGen)
initSpin n g = (take n ( map boolspin list),g)
    where
        list = randoms g :: [Bool]
        boolspin :: Bool -> Int
        boolspin True = 1
        boolspin False = -1

--eager version of iterate, taken from http://stackoverflow.com/questions/8909997/haskell-repeat-a-function-a-large-number-of-times-without-stackoverflow
itereag :: (a -> a) -> a -> [a]
itereag f x = x `seq` (x : itereag f (f x))

main = do
    g <- newStdGen
    let (x0,g') = initSpin 12 g
    let y = [fst (itereag f (x0,g') !! n )| n <- [0..50]]
    let hy = map (ham 0.025 1.00) y 
    mapM_ print $zip y hy