All geek questions in one place

Haskell 1D Random Access Performance

I have a Haskell program that simulates an Ising model with a Metropolis algorithm. The main operation is a screen operation, which takes the sum of the next neighbors in 2D, and then multiplies this by the central element. Then the item may have been updated.

In C ++, where I get decent performance, I use a 1D array and then linearize access to it with simple index arithmetic. In recent months, I took Haskell to expand the horizon, and also tried to implement the Ising model there. The data structure is just a list Bool:

type Spin = Bool
type Lattice = [Spin]

Then I have a fixed degree:

extent = 30

And a function getthat extracts a specific lattice node, including periodic boundary conditions:

-- Wrap a coordinate for periodic boundary conditions.
wrap :: Int -> Int
wrap = flip mod $ extent

-- Converts an unbounded (x,y) index into a linearized index with periodic
-- boundary conditions.
index :: Int -> Int -> Int
index x y = wrap x + wrap y * extent

-- Retrieve a single element from the lattice, automatically performing
-- periodic boundary conditions.
get :: Lattice -> Int -> Int -> Spin
get l x y = l !! index x y

++, , , std::vector .

, get :

COST CENTRE                        MODULE                SRC                       no.     entries  %time %alloc   %time %alloc

         get                       Main                  ising.hs:36:1-26          153     899100    8.3    0.4     9.2    1.9
          index                    Main                  ising.hs:31:1-36          154     899100    0.5    1.2     0.9    1.5
           wrap                    Main                  ising.hs:26:1-24          155          0    0.4    0.4     0.4    0.4
         neighborSum               Main                  ising.hs:(40,1)-(43,56)   133     899100    4.9   16.6    46.6   25.3
          spin                     Main                  ising.hs:(21,1)-(22,17)   135    3596400    0.5    0.4     0.5    0.4
          neighborSum.neighbors    Main                  ising.hs:43:9-56          134     899100    0.9    0.7     0.9    0.7
          neighborSum.retriever    Main                  ising.hs:42:9-40          136     899100    0.4    0.0    40.2    7.6
           neighborSum.retriever.\ Main                  ising.hs:42:32-40         137    3596400    0.2    0.0    39.8    7.6
            get                    Main                  ising.hs:36:1-26          138    3596400   33.7    1.4    39.6    7.6
             index                 Main                  ising.hs:31:1-36          139    3596400    3.1    4.7     5.9    6.1
              wrap                 Main                  ising.hs:26:1-24          141          0    2.7    1.4     2.7    1.4

, Haskell , / , , .

"" , splitAt, ++, , .

- , , ?

:

-- Copyright © 2017 Martin Ueding <dev@martin-ueding.de>

-- Ising model with the Metropolis algorithm. Random choice of lattice site for
-- a spin flip.

import qualified Data.Text
import System.Random

type Spin = Bool
type Lattice = [Spin]

-- Lattice extent is fixed to a square.
extent = 30
volume = extent * extent

temperature :: Double
temperature = 0.0

-- Converts a `Spin` into `+1` or `-1`.
spin :: Spin -> Int
spin True = 1
spin False = (-1)

-- Wrap a coordinate for periodic boundary conditions.
wrap :: Int -> Int
wrap = flip mod $ extent

-- Converts an unbounded (x,y) index into a linearized index with periodic
-- boundary conditions.
index :: Int -> Int -> Int
index x y = wrap x + wrap y * extent

-- Retrieve a single element from the lattice, automatically performing
-- periodic boundary conditions.
get :: Lattice -> Int -> Int -> Spin
get l x y = l !! index x y

-- Computes the sum of neighboring spings.
neighborSum :: Lattice -> Int -> Int -> Int
neighborSum l x y = sum $ map spin $ map retriever neighbors
    where
        retriever = \(x, y) -> get l x y
        neighbors = [(x+1,y), (x-1,y), (x,y+1), (x,y-1)]

-- Computes the energy difference at a certain lattice site if it would be
-- flipped.
energy :: Lattice -> Int -> Int -> Int
energy l x y = 2 * neighborSum l x y * (spin (get l x y))

-- Converts a full lattice into a textual representation.
latticeToString l = unlines lines
    where
        spinToChar :: Spin -> String
        spinToChar True = "#"
        spinToChar False = "."

        line :: String
        line = concat $ map spinToChar l

        lines :: [String]
        lines = map Data.Text.unpack $ Data.Text.chunksOf extent $ Data.Text.pack line

-- Populates a lattice given a random seed.
initLattice :: Int -> (Lattice,StdGen)
initLattice s = (l,rng)
    where
        rng = mkStdGen s

        allRandom :: Lattice
        allRandom = randoms rng

        l = take volume allRandom

-- Performs a single Metropolis update at the given lattice site.
update (l,rng) x y
    | doUpdate = (l',rng')
    | otherwise = (l,rng')
    where
        shift = energy l x y

        r :: Double
        (r,rng') = random rng

        doUpdate :: Bool
        doUpdate = (shift < 0) || (exp (- fromIntegral shift / temperature) > r)

        i = index x y
        (a,b) = splitAt i l
        l' = a ++ [not $ head b] ++ tail b

-- A full sweep through the lattice.
doSweep (l,rng) = doSweep' (l,rng) (extent * extent)

-- Implementation that does the needed number of sweeps at a random lattice
-- site.
doSweep' (l,rng) 0 = (l,rng)
doSweep' (l,rng) i = doSweep' (update (l,rng'') x y) (i - 1)
    where
        x :: Int
        (x,rng') = random rng

        y :: Int
        (y,rng'') = random rng'

-- Creates an IO action that prints the lattice to the screen.
printLattice :: (Lattice,StdGen) -> IO ()
printLattice (l,rng) = do
    putStrLn ""
    putStr $ latticeToString l

dummy :: (Lattice,StdGen) -> IO ()
dummy (l,rng) = do
    putStr "."

-- Creates a random lattice and performs five sweeps.
main = do
    let lrngs = iterate doSweep $ initLattice 2
    mapM_ dummy $ take 1000 lrngs
+4
3

Data.Vector.Unboxed, std::vector. , †. , ST, , , , , Haskell-.

: , , (n) -ish ; . IntMap .

. , Haskell . , Metropolis . , , .

:

neighboursInList :: [a] -> [(a, (Maybe a, Maybe a))]

- map .

-

data Lattice a = Lattice
     { latticeNodes :: [a]
     , latticeLength :: Int }
   deriving (Functor)

data NodeInLattice a = NodeInLattice
     { thisNode :: a
     , xPrev, xNext, yPrev, yNext :: a }
   deriving (Functor)

neighboursInLattice :: Lattice a -> Lattice (NodeInLattice a)

:

  • .
  • .
  • . , repa . , , - , , , , .

† , .

+8

, 5 .

, unboxed vector ( Data.Vector.Unboxed) 1,8 . , System.Random.

random-mersenne-pure64, 0,32 . 0,22 .

, , -, , "" 0,17 .

, unboxed ( , ) , . ( , .)

LCG. .

- extentBits , , ( randomIndex extentBits, , extent).

, True dummy, .

import Data.Bits ((.&.), shiftL)
import Data.Word
import qualified Data.Vector as V

type Spin = Bool
type Lattice = V.Vector Spin

-- Lattice extent is fixed to a square.
extent, extentBits, volume :: Int
extent = 30
extentBits = 5  -- no of bits s.t. 2**5 >= 30
volume = extent * extent

temperature :: Double
temperature = 0.0

-- Converts a `Spin` into `+1` or `-1`.
spin :: Spin -> Int
spin True = 1
spin False = (-1)

-- Wrap a coordinate for periodic boundary conditions.
wrap :: Int -> Int
wrap = flip mod $ extent

-- Converts an unbounded (x,y) index into a linearized index with periodic
-- boundary conditions.
index :: Int -> Int -> Int
index x y = wrap x + wrap y * extent

-- Retrieve a single element from the lattice, automatically performing
-- periodic boundary conditions.
get :: Lattice -> Int -> Int -> Spin
get l x y = l `V.unsafeIndex` index x y

-- Toggle the spin of an element
toggle :: Lattice -> Int -> Int -> Lattice
toggle l x y = l `V.unsafeUpd` [(i, not (l `V.unsafeIndex` i))] -- flip bit at index i
  where i = index x y

-- Computes the sum of neighboring spins.
neighborSum :: Lattice -> Int -> Int -> Int
neighborSum l x y = sum $ map spin $ map (uncurry (get l)) neighbors
    where
        neighbors = [(x+1,y), (x-1,y), (x,y+1), (x,y-1)]

-- Computes the energy difference at a certain lattice site if it would be
-- flipped.
energy :: Lattice -> Int -> Int -> Int
energy l x y = 2 * neighborSum l x y * spin (get l x y)

-- Populates a lattice given a random seed.
initLattice :: Int -> (Lattice,MyGen)
initLattice s = (l, rng')
    where
        rng = newMyGen s
        (allRandom, rng') = go [] rng volume
        go out r 0 = (out, r)
        go out r n = let (a,r') = randBool r
                     in go (a:out) r' (n-1)

        l = V.fromList allRandom

-- Performs a single Metropolis update at the given lattice site.
update :: (Lattice, MyGen) -> Int -> Int -> (Lattice, MyGen)
update (l, rng) x y
  | doUpdate = (toggle l x y, rng')
  | otherwise = (l, rng')
    where
        doUpdate = (shift < 0) || (exp (- fromIntegral shift / temperature) > r)
        shift = energy l x y
        (r, rng') = randDouble rng

-- A full sweep through the lattice.
doSweep :: (Lattice, MyGen) -> (Lattice, MyGen)
doSweep (l, rng) = iterate updateRand (l, rng) !! (extent * extent)

updateRand :: (Lattice, MyGen) -> (Lattice, MyGen)
updateRand (l, rng)
  = let (x, rng') = randIndex rng
        (y, rng'') = randIndex rng'
    in  update (l, rng'') x y

-- Creates a random lattice and performs five sweeps.
main :: IO ()
main = do let lrngs = iterate doSweep (initLattice 2)
              l = fst (lrngs !! 1000)
          print $ V.length (V.filter id l)  -- count the Trues

-- * Random number generation

data MyGen = MyGen Word32

newMyGen :: Int -> MyGen
newMyGen = MyGen . fromIntegral

-- | Get a (positive) integer with given number of bits.
randInt :: Int -> MyGen -> (Int, MyGen)
randInt bits (MyGen s) =
  let s' = 1664525 * s + 1013904223
      mask = (1 `shiftL` bits) - 1
  in  (fromIntegral (s' .&. mask), MyGen s')

-- | Random Bool value
randBool :: MyGen -> (Bool, MyGen)
randBool g = let (i, g') = randInt 1 g
             in  (if i==1 then True else False, g')

-- | Random index
randIndex :: MyGen -> (Int, MyGen)
randIndex g = let (i, g') = randInt extentBits g
              in if i >= extent then randIndex g' else (i, g')

-- | Random [0,1]
randDouble :: MyGen -> (Double, MyGen)
randDouble rng = let (ri, rng') = randInt 32 rng
                 in (fromIntegral ri / (2**32), rng')

MT, , . , randInt, 100%, 100% , , .

import Data.Bits ((.|.), shiftL, shiftR, xor)
import Data.Word
import qualified Data.Vector as V
import System.Random.Mersenne.Pure64

-- replace these definitions:

-- | Mersenne-Twister generator w/ pool of bits
data MyGen = MyGen PureMT !Int !Word64 !Int !Word64

newMyGen :: Int -> MyGen
newMyGen seed = MyGen (pureMT (fromIntegral seed)) 0 0 0 0

-- | Split w into bottom n bits and rest
splitBits :: Int -> Word64 -> (Word64, Word64)
splitBits n w =
  let w2 = w `shiftR` n             -- top 64-n bits
      w1 = (w2 `shiftL` n) `xor` w  -- bottom n bits
  in (w1, w2)

-- | Get a (positive) integer with given number of bits.
randInt :: Int -> MyGen -> (Int, MyGen)
randInt bits (MyGen p lft1 w1 lft2 w2)
  -- generate at least 64 bits
  | let lft = lft1 + lft2, lft < 64
  = let w1' = w1 .|. (w2 `shiftL` lft1)
        (w2', p') = randomWord64 p
    in randInt bits (MyGen p' lft w1' 64 w2')
  | bits > 64 = error "randInt has max of 64 bits"
  -- if not enough bits in first word, get needed bits from second
  | bits > lft1
  = let needed = bits - lft1
        (bts, w2') = splitBits needed w2
        out = (w1 `shiftL` needed) .|. bts
    in (fromIntegral out, MyGen p (lft2 - needed) w2' 0 0)
  -- otherwise, just take enough bits from first word
  | otherwise
  = let (out, w1') = splitBits bits w1
    in (fromIntegral out, MyGen p (lft1 - bits) w1' lft2 w2)
+3

There is another way: a random access list (single or double, or not connected at all). It has access in the worst case (by index) O (log (n)) and does not need ordered data. This is not a Skip List, not an Okasaki List, nor a Haskell List. Performance is amazing (tested up to 2 on 29 elements, half a billion). You can find information about this structure in this article ( https://www.linkedin.com/pulse/random-access-list-armando-giuseppe-bonatto-minella/ )

0
source

Source: https://habr.com/ru/post/1688154/

More articles:


All Articles