arrays - La mejor forma de "hacer un bucle sobre una matriz 2-D", utilizando Repa

performance haskell

Question

Considero que la biblioteca de matrices Repa para Haskell es muy interesante y quería hacer un programa simple para tratar de entender cómo usarlo. También hice una implementación sencilla usando listas, que resultaron ser mucho más rápidas. Mi pregunta principal es cómo podría mejorar el código Repa a continuación para que sea el más eficiente (y espero que también sea muy legible). Soy bastante nuevo usando Haskell, y no pude encontrar ningún tutorial fácilmente comprensible sobre Repa [hay una en Haskell Wiki , que de alguna manera olvidé cuando escribí esto], así que no asumas que sé algo. :) Por ejemplo, no estoy seguro de cuándo usar force o deepSeqArray.

El programa se usa para calcular aproximadamente el volumen de una esfera de la siguiente manera:

1. Se especifican el punto central y el radio de la esfera, así como las coordenadas igualmente espaciadas dentro de un cuboide, que se sabe que abarcan la esfera.
2. El programa toma cada una de las coordenadas, calcula la distancia al centro de la esfera y, si es menor que el radio de la esfera, se usa para sumar el volumen total (aproximado) de la esfera.

A continuación, se muestran dos versiones, una con listas y otra con repa. Sé que el código es ineficiente, especialmente para este caso de uso, pero la idea es hacerlo más complicado más adelante.

Para los siguientes valores, y compilando con "ghc -Odph -fllvm -fforce-recomp -rtsopts -threaded", la versión de la lista tarda 1.4 s, mientras que la versión repa tarda 12 s con + RTS -N1 y 10 s con + RTS - N2, aunque no se han convertido chispas (tengo una máquina Intel de doble núcleo (Core 2 Duo E7400 @ 2.8 GHz) con Windows 7 64, GHC 7.0.2 y llvm 2.8). (Comente la línea correcta en la página principal a continuación para simplemente ejecutar una de las versiones).

¡Gracias por cualquier ayuda!

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

Edición : algunos antecedentes que explican por qué he escrito el código tal como está arriba:

Principalmente escribo código en Matlab, y mi experiencia de mejora del rendimiento proviene principalmente de esa área. En Matlab, normalmente desea realizar sus cálculos utilizando funciones que operan directamente en matrices, para mejorar el rendimiento. Mi implementación del problema anterior, en Matlab R2010b, demora 0.9 segundos con la versión de matriz que se muestra a continuación y 15 segundos con ciclos anidados. Aunque sé que Haskell es muy diferente de Matlab, mi esperanza era que pasar de usar listas a usar arreglos Repa en Haskell mejoraría el rendimiento del código. Las conversiones de listas-> Reparar matrices-> vectores están ahí porque no tengo la habilidad suficiente para reemplazarlas con algo mejor. Por eso pido aportes. :) Por lo tanto, los números de temporización son subjetivos, ya que pueden medir mi rendimiento más que las capacidades de los idiomas, pero ahora es una medida válida, ya que lo que decida qué usaré dependerá de si puedo funciona o no

tl; dr: Entiendo que mi código Repa anterior puede ser estúpido o patológico, pero es lo mejor que puedo hacer en este momento. Me encantaría poder escribir mejor código Haskell, y espero que puedan ayudarme en esa dirección (ya lo hicieron los doctores). :)

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

Edición 2 : versión nueva, más rápida y más simple del código Repa

Ahora he leído un poco más sobre Repa, y he pensado un poco. A continuación se muestra una nueva versión de Repa. En este caso, creo las coordenadas x, y y z como matrices en 3-D, usando la función Repa extender, de una lista de valores (similar a cómo funciona ndgrid en Matlab). Luego mapeo sobre estas matrices para calcular la distancia a la partícula esférica. Finalmente, doblo la matriz de distancia tridimensional resultante, cuento cuántas coordenadas hay dentro de la esfera y luego las multiplico por un factor constante para obtener el volumen aproximado. Mi implementación del algoritmo ahora es mucho más similar a la versión de Matlab anterior, y ya no hay ninguna conversión a vector.

La nueva versión se ejecuta en aproximadamente 5 segundos en mi computadora, una mejora considerable desde arriba. El tiempo es el mismo si uso "subprocesos" al compilar, combinado con "+ RTS -N2" o no, pero la versión de subprocesos supera ambos núcleos de mi computadora. Sin embargo, sí vi algunas gotas de "-N2" ejecutadas a 3,1 segundos, pero no pude reproducirlas más tarde. Tal vez es muy sensible a otros procesos que se ejecutan al mismo tiempo? He cerrado la mayoría de los programas en mi computadora al realizar pruebas comparativas, pero todavía hay algunos programas en ejecución, como los procesos en segundo plano.

Si usamos "-N2" y agregamos el interruptor de tiempo de ejecución para apagar el GC paralelo (-qg), el tiempo se reduce constantemente a ~ 4.1 segundos, y al usar -qa para "usar el sistema operativo para configurar la afinidad del hilo (experimental)", el tiempo se redujo a ~ 3.5 segundos. Mirando el resultado de ejecutar el programa con "+ RTS-s", se realiza mucho menos GC usando -qg.

Esta tarde veré si puedo ejecutar el código en una computadora de 8 núcleos, solo por diversión. :)

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

Perfil del espacio para el nuevo código.

Los mismos tipos de perfiles que Don Stewart hizo a continuación, pero para el nuevo código Repa.

Answer 1

Cambié el código para forzar rcoords y particle_extended , y descubrí que estábamos perdiendo la mayor parte del tiempo dentro de ellos directamente:

COST CENTRE                    MODULE               %time %alloc

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

La mayor mejora individual de este código sería generar esas dos entradas constantes de una mejor manera.

Tenga en cuenta que este es básicamente un algoritmo de transmisión lento, y donde está perdiendo tiempo es el costo irrecuperable de asignar al menos dos arreglos de 24361803 elementos de una sola vez, y luego, probablemente, asignar al menos una o dos veces más o dejar de compartir . El mejor caso para este código, creo, con un optimizador muy bueno y un millón de reglas de reescritura, será emparejar la versión de la lista (que también puede paralelizarse muy fácilmente).

Creo que Dons tiene razón en que Ben & co. Estaré interesado en este punto de referencia, pero mi abrumadora sospecha es que este no es un buen caso de uso para una biblioteca de matrices estricta, y mi sospecha es que Matlab está ocultando algunas optimizaciones inteligentes detrás de su función ngrid (optimizaciones, que concederé, que podría ser útil portar a repa).]

Editar:

Aquí hay una manera rápida y sucia de paralelizar el código de la lista. Importa Control.Parallel.Strategies y luego escribe numberInsideParticles como:

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

Esto muestra una buena aceleración a medida que aumentamos los núcleos (12 segundos en un núcleo a 3.7 segundos en 8), pero la sobrecarga de creación de chispa significa que incluso con 8 núcleos solo coincidimos con la versión no paralela de un solo núcleo. Probé algunas estrategias alternativas y obtuve resultados similares. Una vez más, no estoy seguro de cuánto mejor podemos hacer que una versión de lista de un solo hilo aquí. Dado que los cálculos de cada partícula individual son tan baratos, principalmente estamos enfatizando la asignación, no el cálculo. La gran victoria en algo así, imagino, sería la computación vectorizada más que cualquier otra cosa, y hasta donde sé, eso requiere una codificación manual.

También tenga en cuenta que la versión paralela gasta aproximadamente el 70% de su tiempo en GC, mientras que la versión de un núcleo pasa el 1% de su tiempo allí (es decir, la asignación, en la medida de lo posible, se fusiona eficazmente).

Answer 2

He añadido algunos consejos sobre cómo optimizar los programas Repa en el wiki de Haskell: http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_Tutorial#Optimising_Repa_programs

Answer 3

Notas de revisión de código

• El 47.8% de tu tiempo lo pasas en GC.
• 1.5G se asigna en el montón (!)
• El código de repa parece mucho más complicado que el código de lista.
• Se están produciendo muchos GC paralelos
• Puedo obtener hasta un 300% de eficiencia en una máquina -N4
• Poner más firmas de tipo facilitará el análisis ...
• rsize no se usa (¡parece caro!)
• Usted convierte arrays de repa a vectores, ¿por qué?
• Todos sus usos de (**) podrían ser reemplazados por el más barato (^) en Int .
• Hay un número sospechoso de listas grandes y constantes. Todos tienen que convertirse en matrices, eso parece caro.
• any (==True) es el mismo que or

Perfil de tiempo

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

muestra que su función squared_diff es un poco sospechosa:

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

aunque no veo ninguna solución obvia.

Perfilado espacial

Nada demasiado sorprendente en el perfil del espacio: se ve claramente la fase de lista, luego la fase de vector. La fase de lista asigna mucho, que se recupera.

Al desglosar el montón por tipo, vemos inicialmente que se asignan muchas listas y tuplas (a pedido), luego se asigna y mantiene una gran parte de las matrices:

Una vez más, un poco lo que esperábamos ver ... la matriz no se está asignando especialmente más que el código de la lista (de hecho, un poco menos en general), pero está demorando más en ejecutarse.

Comprobación de fugas de espacio con perfiles de retenedor :

Hay algunas cosas interesantes allí, pero nada sorprendente. zcoords se conserva para la duración de la ejecución del programa de lista, luego algunas matrices (SYSTEM) se asignan para la ejecución de reparación.

Inspeccionar el núcleo

Entonces, en este punto, en primer lugar, asumo que realmente implementó los mismos algoritmos en listas y arreglos (es decir, no se está haciendo ningún trabajo adicional en el caso del arreglo), y no hay una fuga de espacio evidente. Así que mi sospecha es el código de repa mal optimizado. Veamos el núcleo (con ghc-core .

• El código basado en la lista se ve bien.
• El código de matriz parece razonable (es decir, aparecen los primitivos sin caja), pero muy complejo, y mucho de ello.

Alineando todos los CAFs.

Agregué pragmas en línea a todas las definiciones de matriz de nivel superior, con la esperanza de eliminar algunas CAfs, y obtener GHC para optimizar un poco más el código de matriz. Esto realmente hizo que a GHC le costara compilar el módulo (asignando hasta 4.3G y 10 minutos mientras trabajaba en él). Esta es una pista para mí de que GHC no pudo optimizar este programa mucho antes, ya que hay nuevas cosas que hacer cuando aumente los umbrales.

Comportamiento

• El uso de -H puede disminuir el tiempo empleado en GC.
• Intenta eliminar las conversiones de listas a repasas a vectores.
• Todas esas CAF (estructuras de datos constantes de nivel superior) son algo extrañas (un programa real no sería una lista de constantes de nivel superior), de hecho, este módulo es patológico, lo que hace que muchos valores se conserven durante largos períodos. En lugar de ser optimizado lejos. Flotar las definiciones locales hacia adentro.
• Pide ayuda a Ben Lippmeier , el autor de Repa, sobre todo porque hay algunas cosas de optimización funky.