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

Question

我正在尝试使用 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 中的当前范围内指定切片。我考虑过使用 Haskell 中的高效并行模板卷积中讨论的分区，但如果每次迭代都会更改它们，这似乎是一种相当重量级的方法。我还考虑过屏蔽切片（乘以二进制向量） - 但这似乎在大型矩阵上表现很差，因为我将为矩阵中的每个点分配一个屏蔽向量......

我的问题 - 有谁知道是否有计划在 Repa 中添加对范围内切片的支持？或者有没有一种我已经可以解决这个问题的高效方法，也许用不同的方法？

Answer 1

提取子范围是一种索引空间操作，很容易用 fromFunction 表达，尽管我们可能应该为它添加一个更好的包装器到 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]

结果中的元素是通过偏移提供的索引并从源中查找来检索的。这种技术自然地扩展到更高级别的数组。

关于实现并行折叠和扫描，我们将通过在库中添加一个原语来实现。我们不能用 map 来定义并行归约，但我们仍然可以使用延迟数组的整体方法。这将是一个合理的正交扩展。

Answer 2

实际上，我认为主要问题是 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 最近写了 several blog posts 在 Haskell 中导出有效的并行扫描。

积分图像（在二维阵列上）可以通过两次扫描来计算，一次水平扫描，一次垂直扫描。在没有扫描原语的情况下，我实现了一个，效率非常低。

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 原语的有效实现是否可行。