如何在一个范围内使用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中使用高效并行模板卷积中讨论的分区,但如果每次迭代都改变它,这似乎是一个相当重要的方法.我还考虑过屏蔽切片(乘以二进制矢量) - 但这似乎再次对大型矩阵执行得很差,因为我会为矩阵中的每个点分配一个掩码矢量...

我的问题 - 有没有人知道是否有计划增加支持切片范围到维修？或者是否有一种高效的方式我可以解决这个问题,也许采用不同的方法？

Answer 1

提取子范围是一个索引空间操作,很容易用fromFunction表达,尽管我们应该为它添加一个更好的包装器.

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]

通过偏移提供的索引并从源中查找该索引来检索结果中的元素.该技术自然地扩展到更高级别的阵列.

关于实现并行折叠和扫描,我们可以通过向库中添加原语来实现.我们无法在地图方面定义并行缩减,但我们仍然可以使用延迟数组的整体方法.这将是一个合理的正交扩展.

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 最近撰写了几篇 关于在 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 原语进行有效的实现。