arrays - 如何使用Repa在一定范围内获取数组切片

我正在尝试使用Repa来实现累积和功能，以便计算积分图像。我当前的实现如下所示:

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

问题出在elementSlice函数中。在matlab或说numpy中，可以将其指定为array [inner，0:outer]。所以我要寻找的是类似以下内容的东西:

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

但是，似乎不允许在Repa中当前范围内指定 slice 。我考虑过使用Haskell的《高效并行模具卷积》中讨论的分区，但是如果每次迭代都更改分区，这似乎是一种相当重量级的方法。我还考虑过对 slice 进行遮罩(乘以二进制 vector )，但是由于在矩阵中的每个点都将分配遮罩 vector ，因此这似乎在大型矩阵上的表现非常差。

我的问题-有人知道是否有计划增加对Repa范围内 slice 的支持吗？还是有一种我可能已经用一种不同的方法来解决这个问题的高效方法？

最佳答案

实际上，我认为主要问题是Repa没有扫描原语。 (但是，Accelerate库非常相似。)扫描有两种变体，前缀扫描和后缀扫描。给定一维数组
[a_1, ..., a_n]
前缀扫描返回
[0, a_0, a_0 + a_1, ..., a_0 + ... + a_{n-1} ]
而后缀扫描产生
[a_0, a_0 + a_1, ..., a_0 + a_1 + ... + a_n ]
我假设这就是您要使用累计和(cumsum)函数的目的。

前缀和后缀扫描很自然地推广到多维数组，并且具有基于树约简的有效实现。关于该主题的较旧的论文是"Scan Primitives for Vector Computers"。同样，Conal Elliott最近在Haskell中编写了odt_a several blog以获得高效的并行扫描。

可以通过进行两次扫描(水平一次和垂直一次)来计算(在2D阵列上)积分图像。在没有扫描原语的情况下，我的效率很低。

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)

然后可以将用于计算积分图像的函数定义为

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

integralImage = horizScan . vertScan

鉴于已针对Accelerate实施了扫描，将其添加到Repa并不难。我不确定是否可以使用现有的Repa原语进行有效的实现。

关于arrays - 如何使用Repa在一定范围内获取数组切片，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/6170008/