arrays - Cómo tomar una porción de matriz con Repa sobre un rango

haskell

Question

Estoy intentando implementar una función de suma acumulativa usando Repa para calcular imágenes integrales. Mi implementación actual se ve así:

cumsum :: (Elt a, Num a) => Array DIM2 a -> Array DIM2 a
cumsum array = traverse array id cumsum''
where
elementSlice inner outer = slice array (inner :. (0 :: Int))
cumsum'' f (inner :. outer) = Repa.sumAll $ elementSlice inner outer

El problema está en la función elementSlice. En matlab o decir numpy esto podría especificarse como array [inner, 0: outer]. Entonces, lo que estoy buscando es algo como:

slice array (inner :. (Range 0 outer))

Sin embargo, parece que las rodajas no se pueden especificar sobre los rangos actualmente en Repa. Consideré usar particiones como se discutió en Efficient Parallel Stencil Convolution en Haskell, pero esto parece ser un enfoque bastante pesado si se cambiaran cada iteración. También consideré enmascarar el corte (multiplicando por un vector binario), pero de nuevo parecía que funcionaría muy mal en matrices grandes ya que asignaría un vector de máscara para cada punto de la matriz ...

Mi pregunta: ¿alguien sabe si hay planes para agregar soporte para cortar un rango a Repa? ¿O hay una forma efectiva de atacar este problema, tal vez con un enfoque diferente?

Answer 1

En realidad, creo que el principal problema es que Repa no tiene una primitiva de escaneo. (Sin embargo, una biblioteca muy similar a Accelerate sí lo hace). Hay dos variantes de escaneo, escaneo de prefijo y escaneo de sufijo. Dado un conjunto 1D

[a_1, ..., a_n]

regresa el prefijo

[0, a_0, a_0 + a_1, ..., a_0 + ... + a_{n-1} ]

mientras que un escaneo de sufijo produce

[a_0, a_0 + a_1, ..., a_0 + a_1 + ... + a_n ]

Supongo que esto es lo que está buscando con su función acumulativa de suma ( cumsum ).

Los escaneos de prefijo y sufijo se generalizan de manera bastante natural en arreglos multidimensionales y tienen una implementación eficiente basada en la reducción de árboles. Un documento relativamente antiguo sobre el tema es "Primitivas de escaneo para computadoras de Vector" . Además, Conal Elliott ha escrito recientemente varias publicaciones en blogs sobre la obtención de escaneos paralelos eficientes en Haskell.

La imagen integral (en una matriz 2D) puede calcularse haciendo dos escaneos, uno horizontal y otro verticalmente. En ausencia de una primitiva de escaneo, implementé una, muy ineficientemente.

horizScan :: Array DIM2 Int -> Array DIM2 Int
horizScan arr = foldl addIt arr [0 .. n - 1]
where
addIt :: Array DIM2 Int -> Int -> Array DIM2 Int
addIt accum i = accum +^ vs
where
vs = toAdd (i+1) n (slice arr (Z:.All:.i))
(Z:.m:.n) = arrayExtent arr
--
-- Given an @i@ and a length @len@ and a 1D array @arr@
-- (of length n) produces a 2D array of dimensions n X len.
-- The columns are filled with zeroes up to row @i@.
-- Subsequently they are filled with columns equal to the
-- values in @arr.
--
toAdd :: Int -> Int -> Array DIM1 Int -> Array DIM2 Int
toAdd i len arr = traverse arr (/sh -> sh:.(len::Int))
(/_ (Z:.n:.m) -> if m >= i then arr ! (Z:.n) else 0) 

La función para calcular la imagen integral se puede definir como

vertScan :: Array DIM2 Int -> Array DIM2 Int
vertScan = transpose . horizScan . transpose
integralImage = horizScan . vertScan

Dado que el escaneo se ha implementado para Accelerate, no debería ser demasiado difícil agregarlo a Repa. No estoy seguro si es posible o no una implementación eficiente utilizando las primitivas Repa existentes.

Answer 2

La extracción de un subrango es una manipulación del espacio de índice que es bastante fácil de expresar con fromFunction, aunque probablemente deberíamos agregar un envoltorio más agradable para la API.

let arr = fromList (Z :. (5 :: Int)) [1, 2, 3, 4, 5 :: Int]
in fromFunction (Z :. 3) (/(Z :. ix) -> arr ! (Z :. ix + 1))
> [2,3,4]

Los elementos en el resultado se recuperan al compensar el índice proporcionado y buscarlo desde la fuente. Esta técnica se extiende naturalmente a las matrices de mayor rango.

Con respecto a la implementación de pliegues y escaneos paralelos, haríamos esto agregando una primitiva a la biblioteca. No podemos definir reducciones paralelas en términos de mapa, pero aún podemos usar el enfoque general de las matrices demoradas. Sería una extensión razonablemente ortogonal.