Оптимизация производительности числовых массивов в Haskell
Я работаю над алгоритмом генерации ландшафта для мира, подобного MineCraft. В настоящее время я использую симплекс-шум, основанный на реализации в статье "Simplex Noise Demystified" [PDF], поскольку предполагается, что симплекс-шум быстрее и иметь меньше артефактов, чем шум Perlin. Это выглядит довольно прилично (см. Изображение), но пока это также довольно медленно.
![enter image description here]()
Запуск функции шума 10 раз (мне нужен шум с разными длинами волн для таких вещей, как высота поверхности, температура, расположение деревьев и т.д.) с 3 октавами шума для каждого блока в куске (16x16x128 блоков) или около 1 миллиона звонки на шумную функцию в целом, занимает около 700-800 мс. Это, по крайней мере, на порядок медленнее с целью создания местности с любой приличной скоростью, несмотря на то, что в алгоритме нет очевидных дорогостоящих операций (по крайней мере, для меня). Просто пол, по модулю, некоторые поисковые системы и базовая арифметика. Алгоритм (написанный в Haskell) приведен ниже. Комментарии SCC предназначены для профилирования. Я опустил функции 2D-шума, так как они работают одинаково.
g3 :: (Floating a, RealFrac a) => a
g3 = 1/6
{-# INLINE int #-}
int :: (Integral a, Num b) => a -> b
int = fromIntegral
grad3 :: (Floating a, RealFrac a) => V.Vector (a,a,a)
grad3 = V.fromList $ [(1,1,0),(-1, 1,0),(1,-1, 0),(-1,-1, 0),
(1,0,1),(-1, 0,1),(1, 0,-1),(-1, 0,-1),
(0,1,1),( 0,-1,1),(0, 1,-1),( 0,-1,-1)]
{-# INLINE dot3 #-}
dot3 :: Num a => (a, a, a) -> a -> a -> a -> a
dot3 (a,b,c) x y z = a * x + b * y + c * z
{-# INLINE fastFloor #-}
fastFloor :: RealFrac a => a -> Int
fastFloor x = truncate (if x > 0 then x else x - 1)
--Generate a random permutation for use in the noise functions
perm :: Int -> Permutation
perm seed = V.fromList . concat . replicate 2 . shuffle' [0..255] 256 $ mkStdGen seed
--Generate 3D noise between -0.5 and 0.5
simplex3D :: (Floating a, RealFrac a) => Permutation -> a -> a -> a -> a
simplex3D p x y z = {-# SCC "out" #-} 16 * (n gi0 (x0,y0,z0) + n gi1 xyz1 + n gi2 xyz2 + n gi3 xyz3) where
(i,j,k) = {-# SCC "ijk" #-} (s x, s y, s z) where s a = fastFloor (a + (x + y + z) / 3)
(x0,y0,z0) = {-# SCC "x0-z0" #-} (x - int i + t, y - int j + t, z - int k + t) where t = int (i + j + k) * g3
(i1,j1,k1,i2,j2,k2) = {-# SCC "i1-k2" #-} if x0 >= y0
then if y0 >= z0 then (1,0,0,1,1,0) else
if x0 >= z0 then (1,0,0,1,0,1) else (0,0,1,1,0,1)
else if y0 < z0 then (0,0,1,0,1,1) else
if x0 < z0 then (0,1,0,0,1,1) else (0,1,0,1,1,0)
xyz1 = {-# SCC "xyz1" #-} (x0 - int i1 + g3, y0 - int j1 + g3, z0 - int k1 + g3)
xyz2 = {-# SCC "xyz2" #-} (x0 - int i2 + 2*g3, y0 - int j2 + 2*g3, z0 - int k2 + 2*g3)
xyz3 = {-# SCC "xyz3" #-} (x0 - 1 + 3*g3, y0 - 1 + 3*g3, z0 - 1 + 3*g3)
(ii,jj,kk) = {-# SCC "iijjkk" #-} (i .&. 255, j .&. 255, k .&. 255)
gi0 = {-# SCC "gi0" #-} mod (p V.! (ii + p V.! (jj + p V.! kk ))) 12
gi1 = {-# SCC "gi1" #-} mod (p V.! (ii + i1 + p V.! (jj + j1 + p V.! (kk + k1)))) 12
gi2 = {-# SCC "gi2" #-} mod (p V.! (ii + i2 + p V.! (jj + j2 + p V.! (kk + k2)))) 12
gi3 = {-# SCC "gi3" #-} mod (p V.! (ii + 1 + p V.! (jj + 1 + p V.! (kk + 1 )))) 12
{-# INLINE n #-}
n gi (x',y',z') = {-# SCC "n" #-} (\a -> if a < 0 then 0 else
a*a*a*a*dot3 (grad3 V.! gi) x' y' z') $ 0.6 - x'*x' - y'*y' - z'*z'
harmonic :: (Num a, Fractional a) => Int -> (a -> a) -> a
harmonic octaves noise = f octaves / (2 - 1 / int (2 ^ (octaves - 1))) where
f 0 = 0
f o = let r = int $ 2 ^ (o - 1) in noise r / r + f (o - 1)
--Generate harmonic 3D noise between -0.5 and 0.5
harmonicNoise3D :: (RealFrac a, Floating a) => Permutation -> Int -> a -> a -> a -> a -> a
harmonicNoise3D p octaves l x y z = harmonic octaves
(\f -> simplex3D p (x * f / l) (y * f / l) (z * f / l))
Для профилирования я использовал следующий код,
q _ = let p = perm 0 in
sum [harmonicNoise3D p 3 l x y z :: Float | l <- [1..10], y <- [0..127], x <- [0..15], z <- [0..15]]
main = do start <- getCurrentTime
print $ q ()
end <- getCurrentTime
print $ diffUTCTime end start
который дает следующую информацию:
COST CENTRE MODULE %time %alloc
simplex3D Main 18.8 21.0
n Main 18.0 19.6
out Main 10.1 9.2
harmonicNoise3D Main 9.8 4.5
harmonic Main 6.4 5.8
int Main 4.0 2.9
gi3 Main 4.0 3.0
xyz2 Main 3.5 5.9
gi1 Main 3.4 3.4
gi0 Main 3.4 2.7
fastFloor Main 3.2 0.6
xyz1 Main 2.9 5.9
ijk Main 2.7 3.5
gi2 Main 2.7 3.3
xyz3 Main 2.6 4.1
iijjkk Main 1.6 2.5
dot3 Main 1.6 0.7
Чтобы сравнить, я также портировал алгоритм на С#. Производительность там примерно в 3-4 раза быстрее, поэтому я думаю, что я должен делать что-то неправильно. Но даже тогда это не так быстро, как хотелось бы. Поэтому мой вопрос заключается в следующем: может ли кто-нибудь сказать мне, есть ли способы ускорить мою реализацию и/или алгоритм вообще или кто-нибудь знает о другом алгоритме шума, который имеет лучшие характеристики производительности, но похожий вид?
Update:
После выполнения некоторых из предложенных ниже предложений код выглядит следующим образом:
module Noise ( Permutation, perm
, noise3D, simplex3D
) where
import Data.Bits
import qualified Data.Vector.Unboxed as UV
import System.Random
import System.Random.Shuffle
type Permutation = UV.Vector Int
g3 :: Double
g3 = 1/6
{-# INLINE int #-}
int :: Int -> Double
int = fromIntegral
grad3 :: UV.Vector (Double, Double, Double)
grad3 = UV.fromList $ [(1,1,0),(-1, 1,0),(1,-1, 0),(-1,-1, 0),
(1,0,1),(-1, 0,1),(1, 0,-1),(-1, 0,-1),
(0,1,1),( 0,-1,1),(0, 1,-1),( 0,-1,-1)]
{-# INLINE dot3 #-}
dot3 :: (Double, Double, Double) -> Double -> Double -> Double -> Double
dot3 (a,b,c) x y z = a * x + b * y + c * z
{-# INLINE fastFloor #-}
fastFloor :: Double -> Int
fastFloor x = truncate (if x > 0 then x else x - 1)
--Generate a random permutation for use in the noise functions
perm :: Int -> Permutation
perm seed = UV.fromList . concat . replicate 2 . shuffle' [0..255] 256 $ mkStdGen seed
--Generate 3D noise between -0.5 and 0.5
noise3D :: Permutation -> Double -> Double -> Double -> Double
noise3D p x y z = 16 * (n gi0 (x0,y0,z0) + n gi1 xyz1 + n gi2 xyz2 + n gi3 xyz3) where
(i,j,k) = (s x, s y, s z) where s a = fastFloor (a + (x + y + z) / 3)
(x0,y0,z0) = (x - int i + t, y - int j + t, z - int k + t) where t = int (i + j + k) * g3
(i1,j1,k1,i2,j2,k2) = if x0 >= y0
then if y0 >= z0 then (1,0,0,1,1,0) else
if x0 >= z0 then (1,0,0,1,0,1) else (0,0,1,1,0,1)
else if y0 < z0 then (0,0,1,0,1,1) else
if x0 < z0 then (0,1,0,0,1,1) else (0,1,0,1,1,0)
xyz1 = (x0 - int i1 + g3, y0 - int j1 + g3, z0 - int k1 + g3)
xyz2 = (x0 - int i2 + 2*g3, y0 - int j2 + 2*g3, z0 - int k2 + 2*g3)
xyz3 = (x0 - 1 + 3*g3, y0 - 1 + 3*g3, z0 - 1 + 3*g3)
(ii,jj,kk) = (i .&. 255, j .&. 255, k .&. 255)
gi0 = rem (UV.unsafeIndex p (ii + UV.unsafeIndex p (jj + UV.unsafeIndex p kk ))) 12
gi1 = rem (UV.unsafeIndex p (ii + i1 + UV.unsafeIndex p (jj + j1 + UV.unsafeIndex p (kk + k1)))) 12
gi2 = rem (UV.unsafeIndex p (ii + i2 + UV.unsafeIndex p (jj + j2 + UV.unsafeIndex p (kk + k2)))) 12
gi3 = rem (UV.unsafeIndex p (ii + 1 + UV.unsafeIndex p (jj + 1 + UV.unsafeIndex p (kk + 1 )))) 12
{-# INLINE n #-}
n gi (x',y',z') = (\a -> if a < 0 then 0 else
a*a*a*a*dot3 (UV.unsafeIndex grad3 gi) x' y' z') $ 0.6 - x'*x' - y'*y' - z'*z'
harmonic :: Int -> (Double -> Double) -> Double
harmonic octaves noise = f octaves / (2 - 1 / int (2 ^ (octaves - 1))) where
f 0 = 0
f o = let r = 2 ^^ (o - 1) in noise r / r + f (o - 1)
--3D simplex noise
--syntax: simplex3D permutation number_of_octaves wavelength x y z
simplex3D :: Permutation -> Int -> Double -> Double -> Double -> Double -> Double
simplex3D p octaves l x y z = harmonic octaves
(\f -> noise3D p (x * f / l) (y * f / l) (z * f / l))
Вместе с уменьшением размера моего куска до 8x8x128, генерация новых фрагментов рельефа теперь происходит со скоростью около 10-20 кадров в секунду, что означает, что перемещение теперь не так проблематично, как раньше. Конечно, любые другие улучшения производительности по-прежнему приветствуются.
Ответы
Ответ 1
Первоначально это означает, что ваш код очень полиморфен. Вы должны специализировать свой тип с плавающей точкой равномерно на Double, поэтому GHC (и LLVM) имеют возможность применять более агрессивные оптимизации.
Примечание. Для тех, кто пытается воспроизвести, этот код импортирует:
import qualified Data.Vector as V
import Data.Bits
import Data.Time.Clock
import System.Random
import System.Random.Shuffle
type Permutation = V.Vector Int
Ok. Там можно многое сделать, чтобы улучшить этот код.
<сильные > Улучшения
Представление данных
- Специализируется на конкретный тип с плавающей запятой, а не на полиморфные функции с плавающей запятой
- Замените кортеж
(a,a,a) на unboxed triple T !Double !Double !Double
- Переключитесь с
Data.Array на Data.Array.Unboxed на Permutations
- Заменить использование массива с коротким типом с многомерным распакованным массивом из пакета
repa
Флаги компилятора
- Скомпилировать с помощью
-O2 -fvia-C -optc-O3 -fexcess-precision -optc-march=native (или эквивалентно с -fllvm)
- Увеличить порог spec constr -
-fspec-constr-count=16
Более эффективные библиотечные функции
- Используйте mersenne-random вместо StdGen для генерации рандомов.
- Замените
mod на rem
- Замените индексирование
V.! на непроверенную индексацию VU.unsafeIndex (после перехода на Data.Vector.Unboxed
Настройки времени выполнения
- Увеличьте область распределения по умолчанию:
-A20M или -H
Кроме того, проверьте, что ваш алгоритм идентичен С# 1, и вы используете те же структуры данных.
