Я нашел библиотеку массива Repa для Haskell очень интересной и хотел сделать простую программу, чтобы попытаться понять, как ее использовать. Я также сделал простую реализацию, используя списки, которые оказались намного быстрее. Мой главный вопрос: каким образом я могу улучшить код Repa ниже, чтобы сделать его наиболее эффективным (и, надеюсь, также очень читаемым). Я довольно недавно использовал Haskell, и я не мог найти легко понятного учебника по Repa [ edit), который есть в Haskell Wiki, что я как-то забыл, когда писал это], так что не думайте, что я ничего знаю.:) Например, я не уверен, когда использовать силу или deepSeqArray.
Программа используется, чтобы приблизительно рассчитать объем сферы следующим образом:
- Указана центральная точка и радиус сферы, а также одинаково расположенные координаты внутри кубоида, которые, как известно, охватывают сферу.
- Программа берет каждую из координат, вычисляет расстояние до центра сферы и, если она меньше радиуса сферы, она используется для суммирования по общему (приблизительному) объему сферы.
Ниже приведены две версии, одна из которых содержит списки, а другая - с помощью 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 (<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 улучшит производительность кода. Конверсии из списков- > Репа-массивы- > векторы есть, потому что я недостаточно квалифицирован, чтобы заменить их чем-то лучшим. Вот почему я прошу внести свой вклад.:) Таким образом, временные числа, приведенные выше, являются субъективными, поскольку он может измерять мою производительность больше, чем способности языков, но для меня это действительная метрика, так как то, что решает, что я буду использовать, зависит от того, могу ли я сделать это работает или нет.
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 как 3-D массивы, используя функцию расширения Repa, из списка значений (аналогично тому, как ndgrid работает в Matlab). Затем я сопоставляю эти массивы для вычисления расстояния до сферической частицы. Наконец, я сбрасываю полученное трехмерное расстояние, подсчитываю, сколько координат находится внутри сферы, а затем умножьте его на постоянный коэффициент, чтобы получить приблизительный объем. Моя реализация алгоритма теперь намного больше похожа на версию Matlab выше, и больше нет преобразования в вектор.
Новая версия работает примерно через 5 секунд на моем компьютере, что значительно улучшает ее. Сроки совпадают, если я использую "threaded" во время компиляции в сочетании с "+ RTS -N2" или нет, но потоковая версия максимизирует оба ядра моего компьютера. Тем не менее, я видел несколько капель "-N2" до 3,1 секунды, но не смог воспроизвести их позже. Может быть, он очень чувствителен к другим процессам, работающим одновременно? Я закрыл большинство программ на своем компьютере при бенчмарке, но все еще есть некоторые программы, такие как фоновые процессы.
Если мы используем "-N2" и добавляем переключатель времени выполнения для отключения параллельного GC (-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.
