标签: repa

我有一段看起来像这样的Haskell代码:

fst . f $ (Z :. i `div` 2)

Z并:.取自Repa库,定义如下:

data Z = Z deriving (Show, Read, Eq, Ord)
infixl 3 :. 
data tail :. head = !tail :. !head deriving (Show, Read, Eq, Ord)

表达式右侧$定义了一个数组索引,f而是一个获取该索引并返回一对的函数.这将编译为以下Core:

case f_a2pC
       (case ># x_s32E 0 of _ {
          False ->
            case <# x_s32E 0 of _ {
              False -> :. Z (I# (quotInt# x_s32E 2));
              True -> :. Z (I# (-# (quotInt# (+# …

修复和加速API相似性

Haskell repa库用于在CPU上自动并行计算数组.加速库是GPU上的自动数据并行.API非常相似,具有相同的N维数组表示.人们甚至可以切换与加速,维修服务阵列fromRepa和toRepa在Data.Array.Accelerate.IO:

fromRepa :: (Shapes sh sh', Elt e) => Array A sh e -> Array sh' e
toRepa   :: Shapes sh sh'          => Array sh' e  -> Array A sh e

有多个后端用于加速,包括LLVM,CUDA和FPGA(参见http://www.cse.unsw.edu.au/~keller/Papers/acc-cuda.pdf的图2 ).虽然图书馆似乎没有得到维护,但我发现了一个加速的后端.鉴于修复和加速编程模型是相似的,我希望在它们之间有一种优雅的切换方式,即一次写入的函数可以用repa的R.computeP执行,或者用加速的后端执行,例如使用CUDA 运行函数.

两个非常相似的功能:修复和加速南瓜

采用简单的图像处理阈值功能.如果灰度像素值小于50,则将其设置为0,否则保留其值.这是它对南瓜的作用:

以下代码介绍了修复和加速实现:

module Main where

import qualified Data.Array.Repa as R
import qualified Data.Array.Repa.IO.BMP as R
import qualified Data.Array.Accelerate as A
import qualified Data.Array.Accelerate.IO as A
import qualified Data.Array.Accelerate.Interpreter as A

import Data.Word …

我正在使用GHC 7.4.1来尝试编译使用Repa的程序.但是在编译过程中,我的内存耗尽了.有了ghc -v,我可以看到它陷入了SpecConstr阶段.

SpecConstr是GHC的Core-to-Core转换之一.西蒙佩顿琼斯有一个很好的描述在这里,而且有一些代码在这里,但它是相当缓慢的去为我,因为我不是很熟悉,GHC的内部运作.

我希望能够以某种方式帮助编译器 - 有没有办法告诉它哪里被困？或者,有没有办法在这个阶段限制内存使用,直到我可以在更大的机器上重新编译？

谢谢,乍得

更具体地说,我有以下看似无害的小修复3程序:

{-# LANGUAGE QuasiQuotes #-}

import Prelude hiding (map, zipWith)
import System.Environment (getArgs)
import Data.Word (Word8)
import Data.Array.Repa
import Data.Array.Repa.IO.DevIL
import Data.Array.Repa.Stencil
import Data.Array.Repa.Stencil.Dim2

main = do
  [s] <- getArgs
  img <- runIL $ readImage s

  let out = output x where RGB x = img
  runIL . writeImage "out.bmp" . Grey =<< computeP out

output img = map cast . blur . blur $ blur grey
  where
    grey              = traverse img to2D luminance
    cast n            = floor n :: …

我编写了两种算法实现来计算对称矩阵的所有特征值和特征向量.一种实现使用REPA库

https://github.com/felipeZ/Haskell-abinitio/blob/master/Science/QuantumChemistry/NumericalTools/JacobiMethod.hs

而另一个实现使用可变向量和ST monad

https://github.com/felipeZ/Haskell-abinitio/blob/master/Science/QuantumChemistry/NumericalTools/JacobiMethodST.hs

该算法是Jacobi方法,其描述可以在http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm中找到.

我使用100个100 x 100的矩阵测试了这两个实现,顺序运行代码,我发现了以下几次:

                  REPA      Mutable Vectors

Total time(s)    6.7            28.5

GC (s)           0.2             1.2

雅可比算法需要迭代更新一些矩阵条目,这意味着大多数矩阵在迭代之间保持不变.因此,我会猜错(错误地)在REPA实现中为每次迭代复制新矩阵的成本将大于使用monad ST变换矩阵的成本,因为据我所知,REPA不会改变数组但是复制它.

它是REPA融合所有操作并避免在每次迭代中复制新数组？还是别的什么？

有人可以对此结果发表评论吗？

想象一下,我想在一个数组上映射一个函数,但该函数不仅具有类型, a -> b 而且 a -> Int -> b 函数也采用索引.我怎么做？

我正在从磁盘加载RGB图像JuicyPixels-repa.不幸的是,图像的阵列表示是Array F DIM3 Word8内部维度是RGB像素的位置.这与repaRGB图像的现有图像处理算法有点不相容Array U DIM2 (Word8, Word8, Word8).

我想计算图像的RGB直方图,我正在使用签名搜索函数:

type Hist = Array U DIM1 Int
histogram:: Array F DIM3  Word8 -> (Hist, Hist, Hist)

如何折叠我的3d数组以获得每个颜色通道的1d数组？

编辑:

主要的问题不在于我无法从转换DIM3到DIM2每通道(容易切片完成).问题是我必须迭代源图像,DIM2或者DIM3必须累积到DIM1不同Shape (Z:.256) 和范围的数组.所以我不能使用repa,foldS因为它将维度减少了一个,但程度相同.

我也进行了实验,traverse但它迭代了目标图像的范围,提供了从源图像中获取像素的功能,这将导致代码效率非常低,为每个颜色值计算相同的像素.

一个好方法是Vector使用直方图类型作为累加器进行简单的折叠,但不幸的是我没有U(无盒装)或V基于(矢量)的数组,我可以从中有效地得到一个Vector.我有一个Array F(外国指针).

问题

我正在尝试理解Repa如何工作,我正在使用Repa Examples包中的"模糊"示例代码.该代码使用stencil2 Quasi Quote:

[stencil2|   2  4  5  4  2
             4  9 12  9  4
             5 12 15 12  5
             4  9 12  9  4
             2  4  5  4  2 |]

这只是TemplateHaskell片段,它生成一个函数:

makeStencil2 5 5 coeffs where
     {-# INLINE[~0] coeffs #-}
     coeffs = \ ix -> case ix of
                      Z :. -2 :. -2 -> Just 2
                      Z :. -2 :. -1 -> Just 4
                      Z :. -2 :. 0 -> …

在GNU Octave中这段代码 -

[e, ix] = min(X);

将返回最小元素及其位置.你如何修复任意二元函数？

这就是我想出的:

min x = z $ foldl' f (e,0,0) es
  where
    (e:es) = toList x
    f (a,ix,r) b = let ix' = ix+1 in if a < b then (a,ix',r) else (b,ix',ix')
    z (a,ix,r) = (a,r)

在上面的例子中,我们将repa 1D矩阵转换为list并使用foldl'(来自Data.List)和两个累加器 - 一个用于计算迭代(ix),另一个用于保存min元素(r)的位置.但是使用repa的重点是使用数组,而不是列表!

在repa中,Array类型有两个折叠(foldS和foldP) - 但是它们只能取类型函数(a - > a - > a) - 意思是,我不能将带累加器的元组传递给它.还有遍历,原则上可以将1D数组减少为标量数组:

min x = traverse x to0D min
  where
    to0D (Z:.i) = Z
    min f (Z) = ??? -- how to …

最近,我读了纸张在即将到来的广义流的融合vector和DPH图书馆.这似乎是非常有趣的发展.我现在开始尝试DPH(从GHC 7.6开始,并计划在它出来时升级到7.8 SIMD版本).我还从Repa库文档中看到它可以执行并行数组工作.Repa似乎是成熟的版本,相比之下DPH,GHC 7.4没有被认为是准备好的.现在,DPH似乎是成熟的,我想知道什么是之间的主要优点和缺点Repa和DPH封装,GHC 7.6.我搜索了StackOverflow和谷歌,但找不到Repa和之间的比较DPH.因此,这个问题.