Mar*_*ing 4 performance 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.它具有非常快速的随机访问,但它并不真正允许纯功能更新†.你仍然可以通过在STmonad中工作来做这样的更新,事实上这可能是一个能给你带来最佳性能的解决方案,但它并不是真正的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)
这种方法有许多优点:
† 要纯功能更新矢量,您需要制作完整的副本.
| 归档时间： |
|
| 查看次数： |
395 次 |
| 最近记录： |