1D Haskell列表上的随机访问性能

我有一个 Haskell程序,它使用Metropolis模拟Ising模型
算法.主要操作是模板操作,它取下一个的总和
2D中的邻居然后将其与中心元素相乘.那么
元素可能已更新.

在C中,我获得了不错的性能,我使用一维数组然后线性化
使用简单的索引算术访问它.在过去的几个月里,我已经拿起Haskell来拓宽视野,并尝试在那里实施Ising模型.数据结构只是Bool的列表：

type Spin = Bool
type Lattice = [Spin]

然后我有一些固定的程度：

extent = 30

还有一个get函数,用于检索特定的网格站点,包括周期性边界条件：

-- 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

我在C中使用相同的东西,它在那里工作正常,虽然我知道
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

你总是可以使用 Data.Vector.Unboxed,这与std :: vector基本相同.它具有非常快速的随机访问,但它并不真正允许纯功能更新†.你仍然可以通过 ST monad工作来做这样的更新,事实上这可能是给你最佳性能的解决方案,但它并不是真正的Haskell惯用语.

更好：使用一个允许查找和更新以及log(n)-ish时间的功能结构;这是基于树的结构的典型特征. IntMap应该运作得很好.

我不建议这样做.一般来说,在Haskell中,你想要避免任何索引.正如你所说,像Metropolis这样的算法实际上是基于模板.每次旋转的操作都不需要看到比直接邻居更多的操作,因此最好相应地构建程序.

即使在一个简单的列表中,也很容易实现对直接邻居的有效访问：实现

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

那么实际的算法只是这些本地环境的地图.

对于周期性情况,你应该实际上做到这一点

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 library内置了模板支持.在超级计算机上运行的所有代码都必须使用类似的东西,因为访问位于集群中另一个节点上的随机元素比访问处理器自己的节点内存慢得多.

†要纯功能更新矢量,您需要制作完整的副本.