Оптимизация производительности числовых массивов в Haskell - программирование
Подтвердить что ты не робот

Оптимизация производительности числовых массивов в 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 кадров в секунду, что означает, что перемещение теперь не так проблематично, как раньше. Конечно, любые другие улучшения производительности по-прежнему приветствуются.

4b9b3361

Ответ 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, и вы используете те же структуры данных.