Блог

Я нахожу библиотеку массивов Repa для Haskell очень интересной и хотел создать простую программу, чтобы попытаться понять, как ее использовать. Я также создал простую реализацию с использованием списков, которая оказалась намного быстрее. Мой главный вопрос заключается в том, как я мог бы улучшить приведенный ниже код Repa, чтобы сделать его наиболее эффективным (и, надеюсь, также очень читаемым). Я совсем новичок в использовании Haskell, и я не смог найти ни одного легко понятного руководства по Repa [редактировать в Haskell Wiki есть одно, которое я почему-то забыл, когда писал это], так что не думайте, что я что-то знаю. 🙂 Например, я не уверен, когда использовать force или deepSeqArray.

Программа используется для приблизительного вычисления объема сферы следующим образом:

1. Указаны центральная точка и радиус сферы, а также равноудаленные координаты внутри кубоида, которые, как известно, охватывают сферу.
2. Программа берет каждую из координат, вычисляет расстояние до центра сферы, и если оно меньше радиуса сферы, оно используется для суммирования общего (приблизительного) объема сферы.

Ниже показаны две версии, одна с использованием списков, а другая с использованием repa. Я знаю, что код неэффективен, особенно для данного варианта использования, но идея в том, чтобы позже усложнить его.

Для приведенных ниже значений и компиляции с использованием «ghc -Odph -fllvm -fforce-recomp -rtsopts -threaded» для версии списка требуется 1,4 с, в то время как для версии repa требуется 12 с с RTS -N1 и 10 с с RTS -N2, хотя искровые запросы не преобразуются (у меня двухъядерный компьютер Intel (Core 2 Duo E7400 @ 2,8 ГГц) под управлением Windows 7 64, GHC 7.0.2 и llvm 2.8). (Прокомментируйте правильную строку в main ниже, чтобы просто запустить одну из версий.)

Спасибо за любую помощь!

 import Data.Array.Repa as R
import qualified Data.Vector.Unboxed as V
import Prelude as P

-- Calculate the volume of a sphere by putting it in a bath of coordinates. Generate coordinates (x,y,z) in a cuboid. Then, for each coordinate, check if it is inside the sphere. Sum those coordinates and multiply by the coordinate grid step size to find an approximate volume.


particles = [(0,0,0)] -- used for the list alternative --[(0,0,0),(0,2,0)]
particles_repa = [0,0,0::Double] -- used for the repa alternative, can currently just be one coordinate

-- Radius of the particle
a = 4

-- Generate the coordinates. Could this be done more efficiently, and at the same time simple? In Matlab I would use ndgrid.
step = 0.1 --0.05
xrange = [-10,-10 step..10] :: [Double]
yrange = [-10,-10 step..10]
zrange = [-10,-10 step..10]

-- All coordinates as triples. These are used directly in the list version below.
coords = [(x,y,z)  | x <- xrange, y <- yrange, z <- zrange]

---- List code ----

volumeIndividuals = fromIntegral (length particles) * 4*pi*a**3/3

volumeInside = step**3 * fromIntegral (numberInsideParticles particles coords)

numberInsideParticles particles coords = length $ filter (==True) $ P.map (insideParticles particles) coords

insideParticles particles coord =  any (==True) $ P.map (insideParticle coord) particles

insideParticle (xc,yc,zc) (xp,yp,zp) = ((xc-xp)^2 (yc-yp)^2 (zc-zp)^2) < a**2
---- End list code ----

---- Repa code ----

-- Put the coordinates in a Nx3 array.
xcoords = P.map ((x,_,_) -> x) coords
ycoords = P.map ((_,y,_) -> y) coords
zcoords = P.map ((_,_,z) -> z) coords

-- Total number of coordinates
num_coords = (length xcoords) ::Int

xcoords_r = fromList (Z :. num_coords :. (1::Int)) xcoords
ycoords_r = fromList (Z :. num_coords :. (1::Int)) ycoords
zcoords_r = fromList (Z :. num_coords :. (1::Int)) zcoords

rcoords = xcoords_r R.   ycoords_r R.   zcoords_r

-- Put the particle coordinates in an array, then extend (replicate) this array so that its size becomes the same as that of rcoords
particle = fromList (Z :. (1::Int) :. (3::Int)) particles_repa
particle_slice = slice particle (Z :. (0::Int) :. All)
particle_extended = extend (Z :. num_coords :. All) particle_slice

-- Calculate the squared difference between the (x,y,z) coordinates of the particle and the coordinates of the cuboid.
squared_diff = deepSeqArrays [rcoords,particle_extended] ((force2 rcoords) -^ (force2 particle_extended)) **^ 2
(**^) arr pow = R.map (**pow) arr

xslice = slice squared_diff (Z :. All :. (0::Int))
yslice = slice squared_diff (Z :. All :. (1::Int))
zslice = slice squared_diff (Z :. All :. (2::Int))

-- Calculate the distance between each coordinate and the particle center
sum_squared_diff = [xslice,yslice,zslice] `deepSeqArrays` xslice  ^ yslice  ^ zslice

-- Do the rest using vector, since I didn't get the repa variant working.
ssd_vec = toVector sum_squared_diff

-- Determine the number of the coordinates that are within the particle (instead of taking the square root to get the distances above, I compare to the square of the radius here, to improve performance)
total_within = fromIntegral (V.length $ V.filter (<a**2) ssd_vec)
--total_within = foldAll (x acc -> if x < a**2 then acc 1 else acc) 0 sum_squared_diff

-- Finally, calculate an approximation of the volume of the sphere by taking the volume of the cubes with side step, multiplied with the number of coordinates within the sphere.
volumeInside_repa = step**3 * total_within 

-- Helper function that shows the size of a 2-D array.
rsize = reverse . listOfShape . (extent :: Array DIM2 Double -> DIM2)

---- End repa code ----

-- Comment out the list or the repa version if you want to time the calculations separately.
main = do
    putStrLn $ "Step = " P.   show step
    putStrLn $ "Volume of individual particles = " P.   show volumeIndividuals
    putStrLn $ "Volume of cubes inside particles (list) = " P.   show volumeInside
    putStrLn $ "Volume of cubes inside particles (repa) = " P.   show volumeInside_repa
  

Редактировать: Некоторая справочная информация, объясняющая, почему я написал код таким, как указано выше:

В основном я пишу код в Matlab, и мой опыт повышения производительности исходит в основном из этой области. В Matlab вы обычно хотите выполнять свои вычисления, используя функции, работающие непосредственно с матрицами, для повышения производительности. Моя реализация описанной выше проблемы в Matlab R2010b занимает 0,9 секунды при использовании версии матрицы, показанной ниже, и 15 секунд при использовании вложенных циклов. Хотя я знаю, что Haskell сильно отличается от Matlab, я надеялся, что переход от использования списков к использованию Repa-массивов в Haskell улучшит производительность кода. Преобразования из lists-> Repa arrays-> vectors выполняются потому, что я недостаточно опытен, чтобы заменить их чем-то лучшим. Вот почему я запрашиваю ввод. 🙂 Приведенные выше временные показатели, таким образом, субъективны, поскольку они могут измерять мою производительность больше, чем возможности языков, но для меня это верный показатель прямо сейчас, поскольку то, что я буду использовать, зависит от того, смогу ли я заставить это работать или нет.

tl;dr: I understand that my Repa code above may be stupid or pathological, but it’s the best I can do right now. I would love to be able to write better Haskell code, and I hope that you can help me in that direction (dons already did). 🙂

 function archimedes_simple()

particles = [0 0 0]';
a = 4;

step = 0.1;

xrange = [-10:step:10];
yrange = [-10:step:10];
zrange = [-10:step:10];

[X,Y,Z] = ndgrid(xrange,yrange,zrange);
dists2 = bsxfun(@minus,X,particles(1)).^2  ...
    bsxfun(@minus,Y,particles(2)).^2  ...
    bsxfun(@minus,Z,particles(3)).^2;
inside = dists2 < a^2;
num_inside = sum(inside(:));

disp('');
disp(['Step = ' num2str(step)]);
disp(['Volume of individual particles = ' num2str(size(particles,2)*4*pi*a^3/3)]);
disp(['Volume of cubes inside particles = ' num2str(step^3*num_inside)]);

end
  

Edit 2: New, faster and simpler version of the Repa code

Теперь я немного больше прочитал о Repa и немного подумал. Ниже представлена новая версия Repa. В этом случае я создаю координаты x, y и z в виде трехмерных массивов, используя функцию Repa extend, из списка значений (аналогично тому, как ndgrid работает в Matlab). Затем я сопоставляю эти массивы, чтобы вычислить расстояние до сферической частицы. Наконец, я складываю полученный трехмерный массив расстояний, подсчитываю, сколько координат находится внутри сферы, а затем умножаю это на постоянный коэффициент, чтобы получить приблизительный объем. Моя реализация алгоритма теперь намного больше похожа на версию Matlab выше, и больше нет никакого преобразования в вектор.

Новая версия запускается на моем компьютере примерно за 5 секунд, что является значительным улучшением по сравнению с описанным выше. Время остается таким же, если я использую «threaded» при компиляции в сочетании с » RTS -N2″ или нет, но многопоточная версия максимально использует оба ядра моего компьютера. Я, однако, видел несколько сокращений времени выполнения «-N2» до 3,1 секунды, но не смог воспроизвести их позже. Может быть, это очень чувствительно к другим процессам, запущенным одновременно? Я закрыл большинство программ на своем компьютере при тестировании, но некоторые программы все еще запущены, например, фоновые процессы.

Если мы используем «-N2» и добавляем переключатель времени выполнения, чтобы отключить параллельную сборку (-qg), время последовательно сокращается до ~ 4,1 секунды, а с помощью -qa, чтобы «использовать ОС для установки привязки к потоку (экспериментально)», время сокращается до ~ 3,5 секунды. Глядя на результат запуска программы с помощью » RTS -s», гораздо меньше GC выполняется с использованием -qg.

Сегодня днем я посмотрю, смогу ли я запустить код на 8-ядерном компьютере, просто для удовольствия. 🙂

 import Data.Array.Repa as R
import Prelude as P
import qualified Data.List as L

-- Calculate the volume of a spherical particle by putting it in a bath of coordinates.     Generate coordinates (x,y,z) in a cuboid. Then, for each coordinate, check if it is     inside the sphere. Sum those coordinates and multiply by the coordinate grid step size to     find an approximate volume.

particles :: [(Double,Double,Double)]
particles = [(0,0,0)]

-- Radius of the spherical particle
a = 4

volume_individuals = fromIntegral (length particles) * 4*pi*a^3/3

-- Generate the coordinates. 
step = 0.1
coords_list = [-10,-10 step..10] :: [Double]
num_coords = (length coords_list) :: Int

coords :: Array DIM1 Double
coords = fromList (Z :. (num_coords ::Int)) coords_list

coords_slice :: Array DIM1 Double
coords_slice = slice coords (Z :. All)

-- x, y and z are 3-D arrays, where the same index into each array can be used to find a     single coordinate, e.g. (x(i,j,k),y(i,j,k),z(i,j,k)).
x,y,z :: Array DIM3 Double
x = extend (Z :. All :. num_coords :. num_coords) coords_slice
y = extend (Z :. num_coords :. All :. num_coords) coords_slice
z = extend (Z :. num_coords :. num_coords :. All) coords_slice

-- Calculate the squared distance from each coordinate to the center of the spherical     particle.
dist2 :: (Double, Double, Double) -> Array DIM3 Double
dist2 particle = ((R.map (squared_diff xp) x)   (R.map (squared_diff yp) y)   (R.map (    squared_diff zp) z)) 
    where
        (xp,yp,zp) = particle
        squared_diff xi xa = (xa-xi)^2

-- Count how many of the coordinates are within the spherical particle.
num_inside_particle :: (Double,Double,Double) -> Double
num_inside_particle particle = foldAll (acc x -> if x<a^2 then acc 1 else acc) 0 (force     $ dist2 particle)

-- Calculate the approximate volume covered by the spherical particle.
volume_inside :: [Double]
volume_inside = P.map ((*step^3) . num_inside_particle) particles

main = do
    putStrLn $ "Step = " P.   show step
    putStrLn $ "Volume of individual particles = " P.   show volume_individuals
    putStrLn $ "Volume of cubes inside each particle (repa) = " P.   (P.concat . (    L.intersperse ", ") . P.map show) volume_inside

-- As an alternative, y and z could be generated from x, but this was slightly slower in     my tests (~0.4 s).
--y = permute_dims_3D x
--z = permute_dims_3D y

-- Permute the dimensions in a 3-D array, (forward, cyclically)
permute_dims_3D a = backpermute (swap e) swap a
    where
        e = extent a
        swap (Z :. i:. j :. k) = Z :. k :. i :. j
  

Профилирование пространства для нового кода

Те же типы профилей, что и у Дона Стюарта, сделанные ниже, но для нового кода Repa.

Примечания по обзору кода

• 47,8% вашего времени тратится на сборку данных.
• в куче выделено 1,5 Гб (!)
• Код repa выглядит намного сложнее, чем код списка.
• Происходит много параллельного GC
• Я могу добиться эффективности до 300% на машине a -N4
• Добавление большего количества подписей типов упростит анализ…
• rsize не используется (выглядит дорого!)
• Вы преобразуете массивы repa в векторы, зачем?
• Все ваши варианты использования (**) могут быть заменены более дешевым (^) на Int .
• Существует подозрительное количество больших постоянных списков. Все это должно быть преобразовано в массивы — это кажется дорогим.
• any (==True) такой же, как or

Профилирование по времени

 COST CENTRE                    MODULE               %time %alloc

squared_diff                   Main                  25.0   27.3
insideParticle                 Main                  13.8   15.3
sum_squared_diff               Main                   9.8    5.6
rcoords                        Main                   7.4    5.6
particle_extended              Main                   6.8    9.0
particle_slice                 Main                   5.0    7.6
insideParticles                Main                   5.0    4.4
yslice                         Main                   3.6    3.0
xslice                         Main                   3.0    3.0
ssd_vec                        Main                   2.8    2.1
**^                            Main                   2.6    1.4
  

показывает, что ваша функция squared_diff немного подозрительна:

 squared_diff :: Array DIM2 Double
squared_diff = deepSeqArrays [rcoords,particle_extended]
                    ((force2 rcoords) -^ (force2 particle_extended)) **^ 2    
  

хотя я не вижу никакого очевидного исправления.

Профилирование пространства

В профиле пространства нет ничего слишком удивительного: вы четко видите фазу списка, затем векторную фазу. На этапе списка выделяется много ресурсов, которые восстанавливаются.

Разбивая кучу по типам, мы видим, что изначально выделяется множество списков и кортежей (по требованию), затем выделяется и хранится большая часть массивов:

Опять же, примерно то, что мы ожидали увидеть… материал массива выделяет не особенно много, чем код списка (на самом деле, немного меньше в целом), но его выполнение просто занимает намного больше времени.

Проверка на утечки пространства с помощью профилирования фиксатора:

Там есть несколько интересных вещей, но ничего поразительного. zcoords сохраняется на время выполнения программы list, затем некоторые массивы (СИСТЕМНЫЕ) выделяются для выполнения repa.

Проверка ядра

Итак, на данный момент я сначала предполагаю, что вы действительно реализовали одни и те же алгоритмы в списках и массивах (т. Е. В случае массива не выполняется никакой дополнительной работы), и нет очевидной утечки пространства. Итак, я подозреваю, что код repa плохо оптимизирован. Давайте посмотрим на ядро (с ghc-core.

• Код на основе списка выглядит нормально.
• Код массива выглядит разумно (т. Е. появляются распакованные примитивы), но очень сложный, и его много.

Встраивание всех CAFS

Я добавил встроенные программы ко всем определениям массивов верхнего уровня в надежде удалить некоторые из CAFS и заставить GHC немного сложнее оптимизировать код массива. Это действительно заставило GHC с трудом скомпилировать модуль (выделив до 4,3 Гб и 10 минут на работу над ним). Для меня это подсказка, что GHC не смог оптимизировать эту программу задолго до этого, поскольку при увеличении пороговых значений для нее появляются новые возможности.

Действия

• Использование -H может сократить время, затрачиваемое на сборку данных.
• Попробуйте исключить преобразования из списков в repa в векторы.
• Все эти CAFs (постоянные структуры данных верхнего уровня) довольно странные — реальная программа не представляла бы собой список констант верхнего уровня — фактически, этот модуль патологически таков, что множество значений сохраняется в течение длительных периодов времени, вместо того, чтобы быть оптимизированным. Перемещайте локальные определения внутрь.
• Обратитесь за помощью к Бену Липпмайеру, автору Repa, особенно потому, что там происходит какая-то странная оптимизация.

Я изменил код на force rcoords и particle_extended и обнаружил, что мы теряем львиную долю времени непосредственно в них:

 COST CENTRE                    MODULE               %time %alloc

rcoords                        Main                  32.6   34.4
particle_extended              Main                  21.5   27.2
**^                            Main                   9.8   12.7
  

Очевидно, что самым большим улучшением этого кода было бы создание этих двух постоянных входных данных лучшим образом.

Обратите внимание, что это в основном ленивый алгоритм потоковой передачи, и где вы теряете время, так это заниженные затраты на выделение по крайней мере двух массивов из 24361803 элементов за один раз, а затем, вероятно, выделение по крайней мере еще один или два раза или отказ от совместного использования. Я думаю, что самым лучшим вариантом для этого кода с очень хорошим оптимизатором и множеством правил перезаписи будет примерно соответствовать версии списка (которая также может очень легко распараллеливаться).

Я думаю, что донс прав, что Ben amp; co. вас заинтересует этот тест, но я сильно подозреваю, что это не лучший вариант использования библиотеки строгих массивов, и я подозреваю, что matlab скрывает некоторые хитроумные оптимизации за своей ngrid функцией (оптимизации, я согласен, которые было бы полезно перенести в repa).]

Редактировать:

Вот быстрый и грязный способ распараллелить код списка. Импортируйте Control.Parallel.Strategies , а затем запишите numberInsideParticles как:

 numberInsideParticles particles coords = length $ filter id $ 
    withStrategy (parListChunk 2000 rseq) $ P.map (insideParticles particles) coords
  

Это показывает хорошее ускорение по мере увеличения ядер (12 с на одном ядре до 3,7 с на 8), но накладные расходы на создание spark означают, что даже для 8 ядер мы подбираем только одноядерную непараллельную версию. Я попробовал несколько альтернативных стратегий и получил похожие результаты. Опять же, я не уверен, насколько лучше мы можем сделать, чем однопоточная версия списка здесь. Поскольку вычисления для каждой отдельной частицы очень дешевы, мы в основном делаем упор на распределение, а не на вычисления. Большой выигрыш в чем-то подобном, я полагаю, был бы в векторизованных вычислениях больше, чем в чем-либо другом, и, насколько я знаю, это в значительной степени требует ручного кодирования.

Также обратите внимание, что параллельная версия тратит примерно 70% своего времени на GC, в то время как одноядерная версия тратит там 1% своего времени (т. Е. распределение, насколько это возможно, эффективно распределяется.).

Я добавил несколько советов о том, как оптимизировать программы Repa в Haskell wiki:http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_Tutorial#Optimising_Repa_programs