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

Question

我有一个 Haskell 程序，它用 Metropolis 模拟 Ising 模型 算法。主要操作是一个模板操作，它取 next 的总和 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 这样的算法实际上是基于 stencil。每次自旋的操作应该只需要查看其直接邻居，因此最好相应地构建您的程序。

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

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

---

^†_{要纯功能更新向量，您需要制作完整的副本。}

Answer 2

关闭分析后，您的原始版本在我的笔记本电脑上运行大约需要 5 秒。

将代码转换为使用不可变的、未装箱的向量（来自Data.Vector.Unboxed）是一个简单的修改，并将运行时间减少到大约 1.8 秒。分析该版本表明时间由非常慢的 System.Random 生成器控制。

使用基于 random-mersenne-pure64 包的自定义生成器，我可以将运行时间缩短到大约 0.32 秒。使用线性同余生成器可以将时间缩短到 0.22 秒。

重新分析，瓶颈似乎是对向量操作的边界检查，因此将它们替换为“不安全”的对应物可以将运行时间缩短到大约 0.17 秒。

此时，转换为可变的、未装箱的向量（这是一个比以前更复杂的修改）并没有明显提高性能，但我并没有努力优化它。 （我见过其他算法从使用可变向量中受益匪浅。）

我的 LCG 版本的最终代码如下。我试图在合理的范围内保留尽可能多的原始代码。

一个烦人的地方是必须为随机索引生成指定extentBits，请注意，如果范围是 2 的幂，算法将是最有效的（因为randomIndex 使用给定的数字生成索引extentBits，然后重试直到索引小于extent）。

请注意，我决定在最终格中打印 Trues 的数量，而不是使用 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)