I recently learned Haskell, and I am trying to apply it to the code I use in order to get a feeling for the language.
I really like the Repa library since I manipulate a lot of multi-dimensional data. Yet it is, in my sense, missing a lot of signal processing tools. For instance, I use quite a lot of Gaussian filter combinations, and approximate them with the Canny-Deriche recursive filter, so that the complexity remains linear and independent of the deviation.
So, this is the code that I have started to write, only for Gaussian convolutions, but in any dimensions.
import Control.Monad
import Data.List (tails)
import Data.Array.Repa as R hiding ((++), map, zipWith)
import Data.Array.Repa.Eval
import Data.Array.Repa.Eval.Gang
import Data.Array.Repa.Repr.Unboxed
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
-- |Apply a recursive Canny-Deriche filter on one dimension of a Repa array.
deriche :: Shape sh =>
Double -> Int -> Int -> Bool ->
Array U sh Double ->
IO (Array U sh Double)
deriche s o d bc a
-- Convolution with a dirac is a nop
| s == 0 = return a
-- Convolution only make sense for positive sigmas
| s < 0 = error "Sigma for the filter must be positive"
-- Make sure the selected dimension exists
| d >= rank (extent a) =
error "Invalid array dimension selected for the filter"
| otherwise =
case o of
0 -> let k = (1-ena) * (1-ena) / (1 + 2*alpha*ena - ens)
a0 = k
a1 = k * ena * (alpha - 1)
a2 = k * ena * (alpha + 1)
a3 = -k * ens
par = 1
in do
b <- UM.new $ size $ extent a :: IO (UM.IOVector Double)
innerLoop a0 a1 a2 a3 0 [] True b
liftM (fromUnboxed $ extent a) (U.unsafeFreeze b)
otherwise -> error "Unimplemented filter order"
where alpha = 1.695 / s
expAlpha = exp alpha
ena = exp (-alpha)
ens = ena * ena
b1 = -2 * ena
b2 = ens
ndim = rank sha
dims = reverse $ listOfShape sha
ddim = dims!!d
stride = product $ tails dims !! (d+1)
sha = extent a
innerLoop :: Double -> Double -> Double -> Double ->
Int -> [Int] -> Bool -> UM.IOVector Double ->
IO ()
innerLoop a0 a1 a2 a3 l path pl b
| l < ndim && l == d = innerLoop a0 a1 a2 a3 (l+1) (0:path) pl b
| l < ndim && pl =
let len = dims!!l
threads = gangSize theGang
chunkLen = len `quot` threads
chunkOver = len `rem` threads
splitIdx thread
| thread < chunkOver = thread * (chunkLen + 1)
| otherwise = thread * chunkLen + chunkOver
fill :: Int -> Int -> IO ()
fill ix end =
mapM_ (\i -> innerLoop a0 a1 a2 a3 (l+1) (i:path) False b) [ix..end-1]
in gangIO theGang $ \thread ->
let start = splitIdx thread
end = splitIdx (thread + 1)
in fill start end
| l < ndim =
mapM_ (\i -> innerLoop a0 a1 a2 a3 (l+1) (i:path) pl b) [0..dims!!l-1]
| otherwise = do
-- First pass
if bc
then undefined --pass1 0 (a `linearIndex` idx)
else pass1 0 idx 0 0 0
-- Second pass
if bc
then undefined -- pass2 ...
else pass2 (ddim-1) (idx+(ddim-1)*stride) 0 0 0 0
where idx = toIndex sha (shapeOfList $ cop path)
cop = id
pass1 !i !j !xp !yp !yb
| i == ddim = return ()
| otherwise = do
UM.unsafeWrite b j yc
pass1 (i+1) (j+stride) xc yc yp
where yc = a0*xc + a1*xp - b1*yp - b2*yb
xc = a `unsafeLinearIndex` j
{-# INLINE pass1 #-}
pass2 !i !j !xn !xa !yn !ya
| i < 0 = return ()
| otherwise = do
bj <- UM.unsafeRead b j
UM.unsafeWrite b j (bj + yc)
pass2 (i-1) (j-stride) xc xn yc yn
where yc = a2*xn + a3*xa - b1*yn - b2*ya
xc = a `unsafeLinearIndex` j
{-# INLINE pass2 #-}
{-# INLINE innerLoop #-}
-- |Sequentially apply a Canny-Deriche filter on two dimensions.
deriche2D :: Double -> Int -> Bool ->
Array U DIM2 Double ->
IO (Array U DIM2 Double)
deriche2D s o bc = deriche s o 0 bc >=> deriche s o 1 bc
Since I am a complete novice in Haskell, I would appreciate any comment you might give. In particular the code I developed is quite un-repa-ish due to the recursive filter. Also I have no clue as to when I should use bangs, since adding them to any array variable has no positive effect. Any idea on how to share the arrays for deriche2D
would be greatly appreciated.
Also, I am compiling with -rtsopts -O2 -threaded -fllvm
and running with +RTS -N4
.