Лучший способ «зацикливаться на двумерном массиве», используя Repa

Я нахожу библиотеку массивов Repa для Haskell очень интересной, и я хотел создать простую программу, чтобы попытаться понять, как ее использовать. Я также сделал простую реализацию с использованием списков, которая оказалась намного быстрее. Мой главный вопрос заключается в том, как я могу улучшить приведенный ниже код Repa, чтобы сделать его наиболее эффективным (и, надеюсь, также очень читаемым). Я новичок в использовании Haskell, и я не мог Я не нашел ни одного легко понятного руководства по Repa [ edit , есть одно в Haskell Wiki , которое я как-то забыл, когда писал это], поэтому не думайте, что я что-то знаю. :) Например, я не уверен, когда использовать силу или 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). (Закомментируйте правильную строку в основном тексте ниже, чтобы просто запустить одну из версий.)

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

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 ( 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 улучшит производительность кода. Преобразования из списков-> массивов Repa-> векторов существуют, потому что я недостаточно опытен, чтобы заменить их чем-то лучшим. Вот почему я прошу вас внести свой вклад. :) Приведенные выше временные числа, таким образом, являются субъективными, поскольку они могут измерять мою производительность больше, чем способности языков, но это действительный показатель для меня прямо сейчас, поскольку от того, что решает, что я буду использовать зависит от того, могу ли я заставить его работать или нет.

tl; dr: Я понимаю, что приведенный выше код Repa может быть глупым или патологическим, но это лучшее, что я могу сделать прямо сейчас. Я хотел бы иметь возможность писать лучший код на Haskell, и я надеюсь, что вы можете помочь мне в этом направлении (доны уже сделали). :)

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

Редактировать 2 : Новая, более быстрая и простая версия кода Repa

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

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

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

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

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