haskell - Mal desempeño con transposición y suma acumulada en Repa.

Question

He desarrollado una función de suma acumulativa como se define a continuación en la biblioteca de Haskell Repa. Sin embargo, me he encontrado con un problema al combinar esta función con la operación de transposición. Las 3 operaciones siguientes toman menos de un segundo:

cumsum $ cumsum $ cumsum x
transpose $ transpose $ transpose x
transpose $ cumsum x

Sin embargo, si escribo:

cumsum $ transpose x

El rendimiento se degrada horriblemente. Mientras que cada operación individual en aislamiento toma menos de un segundo en una imagen de 1920x1080, cuando se combinan ahora toman más de 30 segundos ...

¿Alguna idea de qué puede estar causando esto? Mi instinto me dice que tiene algo que ver con arreglos retrasados, no forzar en el momento adecuado, etc. Pero todavía no tengo suficiente experiencia para rastrear esto.

{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}

import Data.Array.Repa as Repa

{-# INLINE indexSlice #-}
indexSlice :: (Shape sh, Elt a) => Int -> Array (sh :. Int) a -> (sh :. Int) -> a
indexSlice from arr (z :. ix) = arr `unsafeIndex` (z :. (ix + from))

{-# INLINE sliceRange #-}
sliceRange :: (Slice sh, Shape sh, Elt a) => Int -> Int -> Array (sh :. Int) a -> Array (sh :. Int) a
sliceRange from to arr = fromFunction (z :. (to - from + 1)) $ indexSlice from arr
    where (z :. _) = extent arr

{-# INLINE cumsum'' #-}
cumsum'' :: (Slice (SliceShape sh), Slice sh, Shape (FullShape sh), Shape (SliceShape sh), Elt a, Num a) =>
     Array (FullShape sh :. Int) a -> t -> (sh :. Int) -> a
cumsum'' arr f (sh :. outer) = Repa.sumAll $ sliceRange 0 outer $ Repa.slice arr (sh :. All)

{-# INLINE cumsum #-}
cumsum :: (FullShape sh ~ sh, Slice sh, Slice (SliceShape sh), Shape sh, Shape (SliceShape sh), Elt a, Num a) =>
    Array (sh :. Int) a -> Array (sh :. Int) a
cumsum arr = Repa.force $ unsafeTraverse arr id $ cumsum'' arr

Answer 1

Desde la perspectiva de un implementador de bibliotecas, la forma de depurar esto es crear un envoltorio para la operación sospechosa, luego mirar el código central para ver si la fusión ha funcionado.

-- Main.hs ---------------------------------------------------
import Solver
import Data.Array.Repa.IO.BMP

main 
 = do   Right img       <- readImageFromBMP "whatever.bmp"
        print $ cumsumBMP img

-- Solver.hs --------------------------------------------------
{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}
module Solver (cumsumBMP) where
import Data.Array.Repa  as Repa
import Data.Word

{- all your defs -}

{-# NOINLINE cumsumBMP #-}
cumsumBMP :: Array DIM3 Word8 -> Array DIM3 Word8
cumsumBMP img = cumsum $ transpose img

He puesto el código "solucionador" en un módulo separado, por lo que solo tenemos que revisar el código central para las definiciones que nos interesan.

Compilar como

touch Solver.hs ; ghc -O2 --make Main.hs /
 -ddump-simpl -dsuppress-module-prefixes -dsuppress-coercions  > dump

Vaya a la definición de cumsumBMP y busque la palabra clave letrec . La búsqueda de letrec es una forma rápida de encontrar los bucles internos.

No muy abajo veo esto: (reformateo ligeramente)

case gen_a1tr
of _ {
  GenManifest vec_a1tv ->
    case sh2_a1tc  `cast` ... of _ { :. sh3_a1iu  sh4_a1iv ->
    case ix''_a1t9  `cast` ... of _ { :. sh1''_a1iz sh2''_a1iA ->
    case sh3_a1iu  `cast` ... of _ { :. sh5_X1n0  sh6_X1n2 ->
    case sh1''_a1iz `cast` ... of _ { :. sh1''1_X1n9 sh2''1_X1nb ->
    case sh5_X1n0             of _ { :. sh7_X1n8   sh8_X1na ->
    ...
    case sh2''1_X1nb           of _ { I# y3_X1nO ->
    case sh4_a1iv             of _ { I# y4_X1nP ->
    case sh2''_a1iA            of _ { I# y5_X1nX ->
    ...
    let { x3_a1x6 :: Int# [LclId]
      x3_a1x6 =
        +#
          (*#
             (+#
                (*#
                   y1_a1iM
                   y2_X1nG)
                y3_X1nO)
             y4_X1nP)
          y5_X1nX } in
    case >=#
           x3_a1x6
           0
    of ...

¡Desastre! El enlace x3_a1x6 está haciendo claramente un trabajo útil (multiplicaciones, adiciones y similares) pero está envuelto en una larga serie de operaciones de desempaquetado que también se ejecutan para cada iteración de bucle. Lo que es peor es que es unboxing la longitud y el ancho (forma) de la matriz en cada iteración, y esta información siempre será la misma. GHC debería realmente flotar estas expresiones de casos fuera del bucle, pero aún no lo hace. Esta es una instancia del problema nº 4081 en el trac de GHC , que se espera que se solucione pronto.

La deepSeqArray es aplicar deepSeqArray a la matriz entrante. Esto coloca una demanda en su valor en el nivel superior (fuera del bucle), lo que le permite a GHC saber que está bien mover las coincidencias de casos más hacia arriba. Para una función como cumsumBMP , también esperamos que la matriz entrante ya se manifieste, por lo que podemos agregar una coincidencia de caso explícita para esto:

{-# NOINLINE cumsumBMP #-}
cumsumBMP :: Array DIM3 Word8 -> Array DIM3 Word8
cumsumBMP img@(Array _ [Region RangeAll (GenManifest _)])
  = img `deepSeqArray` cumsum $ transpose img

Compilando de nuevo, el bucle interno ahora se ve mucho mejor:

letrec {
$s$wfoldlM''_loop_s2mW [...]
  :: Int# -> Word# -> Word# [...]
$s$wfoldlM''_loop_s2mW =
  / (sc_s2mA :: Int#) (sc1_s2mB :: Word#) ->
    case <=# sc_s2mA a_s2ji of _ {
      False -> sc1_s2mB;
      True ->
        $s$wfoldlM''_loop_s2mW
          (+# sc_s2mA 1)
          (narrow8Word#
             (plusWord#
                sc1_s2mB
                (indexWord8Array#
                   rb3_a2gZ
                   (+#
                      rb1_a2gX
                      (+#
                         (*#
                            (+#
                               (*#
                                  wild19_X1zO
                                  ipv1_X1m5)
                               sc_s2mA)
                            ipv2_X1m0)
                         wild20_X1Ct)))))
    }; } in

Eso es un bucle recursivo cerrado y de cola que solo usa operaciones primitivas. Siempre que compile con -fllvm -optlo-O3 , no hay razón para que no se ejecute tan rápido como un programa de C equivalente.

Sin embargo, hay un pequeño hipo al ejecutarlo:

desire:tmp benl$ ./Main 
Main: Solver.hs:(50,1)-(51,45): Non-exhaustive patterns in function cumsumBMP

Esto nos recuerda que debemos forzar la matriz antes de llamar a cumsumBMP .

-- Main.hs ---------------------------------------------------
...
import Data.Array.Repa as Repa
main 
 = do   Right img       <- readImageFromBMP "whatever.bmp"
        print $ cumsumBMP $ Repa.force img

En resumen:

1. Es necesario agregar algunos elementos de deepSeqArray y de coincidencia de patrones a sus funciones de nivel superior para deepSeqArray una deepSeqArray actual en GHC. Esto se demuestra en la versión final de la función cumsumBMP anterior. Si desea que GHC HQ arregle esto pronto, entonces agréguese a sí mismo como un cc al Problema n.º 4081 en el tracto de GHC . Los programas Repa serán mucho más bonitos cuando esto se arregle.
2. No es necesario agregar el pegote a cada función. En este ejemplo no tuve que tocar a indexSlice y amigos. La regla general es agregar el goop a las funciones que usan force , fold o sumAll . Estas funciones instancian los bucles reales que operan sobre los datos de la matriz, es decir, convierten una matriz retrasada en un valor manifiesto.
3. El rendimiento de una parte del código Repa está determinado tanto por el contexto en el que se utiliza como el código real. Si pasa sus arrays retrasados ​​de funciones de nivel superior, entonces se ejecutarán muy lentamente. Hay más discusión sobre esto en el tutorial de Repa .
4. Los archivos BMP que se leen con la biblioteca de repa-io no están pre-forzados, por lo que debe forzarlos antes de usarlos. Este es probablemente el valor predeterminado incorrecto, por lo que lo cambiaré en la próxima versión.