performance - 一维 Haskell 列表上的随机访问性能

Question

我有一个 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

Answer 1

您始终可以使用 Data.Vector.Unboxed ，它与 std::vector 基本相同。它具有非常快的随机访问，但它并不真正允许纯功能更新†。您仍然可以通过在 ST monad 中工作来进行此类更新，实际上这可能是为您提供最佳性能的解决方案，但它并不是真正的 Haskell 惯用语。

更好:使用允许查找和更新以及 log(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 library 内置了模板支持。所有在 super 计算机上运行的代码都必须使用类似的东西，因为访问位于集群中另一个节点上的随机元素比访问处理器自己的节点内存要慢得多。

†要以纯功能方式更新向量，您需要制作完整的副本。