diff --git a/massiv.cabal b/massiv.cabal --- a/massiv.cabal +++ b/massiv.cabal @@ -1,5 +1,5 @@ name: massiv -version: 0.2.7.0 +version: 0.2.8.0 synopsis: Massiv (Массив) is an Array Library. description: Multi-dimensional Arrays with fusion, stencils and parallel computation. homepage: https://github.com/lehins/massiv diff --git a/src/Data/Massiv/Array/Stencil.hs b/src/Data/Massiv/Array/Stencil.hs --- a/src/Data/Massiv/Array/Stencil.hs +++ b/src/Data/Massiv/Array/Stencil.hs @@ -49,10 +49,11 @@ where !window = Window - { windowStart = sCenter - , windowSize = liftIndex2 (-) sz (liftIndex2 (-) sSz (pureIndex 1)) + { windowStart = liftIndex2 min sCenter (liftIndex (max 0) (liftIndex (subtract 1) sz)) + , windowSize = liftIndex (max 0) (liftIndex2 min windowSz (liftIndex2 (-) sz sCenter)) , windowIndex = unValue . stencilF (Value . unsafeIndex arr) } + !windowSz = liftIndex (max 0) (liftIndex2 (-) sz (liftIndex (subtract 1) sSz)) !sz = size arr {-# INLINE mapStencil #-} diff --git a/src/Data/Massiv/Array/Stencil/Convolution.hs b/src/Data/Massiv/Array/Stencil/Convolution.hs --- a/src/Data/Massiv/Array/Stencil/Convolution.hs +++ b/src/Data/Massiv/Array/Stencil/Convolution.hs @@ -2,7 +2,7 @@ {-# LANGUAGE BangPatterns #-} -- | -- Module : Data.Massiv.Array.Stencil.Convolution --- Copyright : (c) Alexey Kuleshevich 2018 +-- Copyright : (c) Alexey Kuleshevich 2018-2019 -- License : BSD3 -- Maintainer : Alexey Kuleshevich -- Stability : experimental @@ -28,20 +28,23 @@ -- Here is how to create a 2D horizontal Sobel Stencil: -- -- > sobelX :: Num e => Stencil Ix2 e e --- > sobelX = makeConvolutionStencil (3 :. 3) (1 :. 1) $ --- > \f -> f (-1 :. -1) 1 . f (-1 :. 1) (-1) . --- > f ( 0 :. -1) 2 . f ( 0 :. 1) (-2) . --- > f ( 1 :. -1) 1 . f ( 1 :. 1) (-1) +-- > sobelX = makeConvolutionStencil (Sz2 3 3) (1 :. 1) $ +-- > \f -> f (-1 :. -1) (-1) . f (-1 :. 1) 1 . +-- > f ( 0 :. -1) (-2) . f ( 0 :. 1) 2 . +-- > f ( 1 :. -1) (-1) . f ( 1 :. 1) 1 -- > {-# INLINE sobelX #-} -- +-- @since 0.1.0 makeConvolutionStencil :: (Index ix, Num e) => Sz ix -> ix -> ((ix -> Value e -> Value e -> Value e) -> Value e -> Value e) -> Stencil ix e e -makeConvolutionStencil !sSz !sCenter relStencil = validateStencil 0 $ Stencil sSz sCenter stencil +makeConvolutionStencil !sz !sCenter relStencil = + validateStencil 0 $ Stencil sz sInvertCenter stencil where + !sInvertCenter = liftIndex2 (-) (liftIndex (subtract 1) sz) sCenter stencil getVal !ix = (inline relStencil $ \ !ixD !kVal !acc -> getVal (liftIndex2 (-) ix ixD) * kVal + acc) 0 {-# INLINE stencil #-} @@ -51,24 +54,30 @@ -- | Make a stencil out of a Kernel Array. This `Stencil` will be slower than if -- `makeConvolutionStencil` is used, but sometimes we just really don't know the -- kernel at compile time. +-- +-- @since 0.1.0 makeConvolutionStencilFromKernel :: (Manifest r ix e, Num e) => Array r ix e -> Stencil ix e e -makeConvolutionStencilFromKernel kArr = Stencil sz sCenter stencil +makeConvolutionStencilFromKernel kArr = Stencil sz sInvertCenter stencil where - !sz = size kArr - !sCenter = liftIndex (`div` 2) sz + !sz@(Sz szi) = size kArr + !szi1 = liftIndex (subtract 1) szi + !sInvertCenter = liftIndex2 (-) szi1 sCenter + !sCenter = liftIndex (`quot` 2) szi stencil getVal !ix = Value (ifoldlS accum 0 kArr) where + !ixOff = liftIndex2 (+) ix sCenter accum !acc !kIx !kVal = - unValue (getVal (liftIndex2 (+) ix (liftIndex2 (-) sCenter kIx))) * kVal + acc + unValue (getVal (liftIndex2 (-) ixOff kIx)) * kVal + acc {-# INLINE accum #-} {-# INLINE stencil #-} {-# INLINE makeConvolutionStencilFromKernel #-} - --- | Make a stencil. +-- | Make a stencil +-- +-- @since 0.1.5 makeCorrelationStencil :: (Index ix, Num e) => Sz ix @@ -82,10 +91,11 @@ {-# INLINE stencil #-} {-# INLINE makeCorrelationStencil #-} - --- | Make a stencil out of a Kernel Array. This `Stencil` will be slower than if --- `makeCorrelationStencil` is used, but sometimes we just really don't know the --- kernel at compile time. +-- | Make a stencil out of a +-- Kernel Array. This `Stencil` will be slower than if `makeCorrelationStencil` is used, but +-- sometimes we just really don't know the kernel at compile time. +-- +-- @since 0.1.5 makeCorrelationStencilFromKernel :: (Manifest r ix e, Num e) => Array r ix e @@ -95,8 +105,9 @@ !sz = size kArr !sCenter = liftIndex (`div` 2) sz stencil getVal !ix = Value (ifoldlS accum 0 kArr) where + !ixOff = liftIndex2 (-) ix sCenter accum !acc !kIx !kVal = - unValue (getVal (liftIndex2 (+) ix (liftIndex2 (+) sCenter kIx))) * kVal + acc + unValue (getVal (liftIndex2 (+) ixOff kIx)) * kVal + acc {-# INLINE accum #-} {-# INLINE stencil #-} {-# INLINE makeCorrelationStencilFromKernel #-} diff --git a/src/Data/Massiv/Array/Stencil/Unsafe.hs b/src/Data/Massiv/Array/Stencil/Unsafe.hs --- a/src/Data/Massiv/Array/Stencil/Unsafe.hs +++ b/src/Data/Massiv/Array/Stencil/Unsafe.hs @@ -48,10 +48,11 @@ where !window = Window - { windowStart = sCenter - , windowSize = liftIndex2 (-) sz (liftIndex2 (-) sSz (pureIndex 1)) + { windowStart = liftIndex2 min sCenter (liftIndex (max 0) (liftIndex (subtract 1) sz)) + , windowSize = liftIndex (max 0) (liftIndex2 min windowSz (liftIndex2 (-) sz sCenter)) , windowIndex = stencil (Just . unsafeIndex arr) } + !windowSz = liftIndex (max 0) (liftIndex2 (-) sz (liftIndex (subtract 1) sSz)) stencil getVal !ix = inline relStencil $ \ !ixD -> getVal (liftIndex2 (+) ix ixD) {-# INLINE stencil #-} !sz = size arr diff --git a/tests/Data/Massiv/Array/StencilSpec.hs b/tests/Data/Massiv/Array/StencilSpec.hs --- a/tests/Data/Massiv/Array/StencilSpec.hs +++ b/tests/Data/Massiv/Array/StencilSpec.hs @@ -1,18 +1,19 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MonoLocalBinds #-} -{-# LANGUAGE OverloadedLists #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MonoLocalBinds #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE OverloadedLists #-} +{-# LANGUAGE ScopedTypeVariables #-} module Data.Massiv.Array.StencilSpec (spec) where -import Control.DeepSeq (deepseq) -import Data.Massiv.CoreArbitrary as A -import Data.Maybe (fromJust) -import Data.Proxy -import Test.Hspec -import Test.QuickCheck -import Test.QuickCheck.Function -import Data.Default (Default(def)) +import Control.DeepSeq (deepseq) +import Data.Default (Default(def)) +import Data.Massiv.CoreArbitrary as A +import Data.Maybe (fromJust) +import Data.Proxy +import Test.Hspec +import Test.QuickCheck +import Test.QuickCheck.Function -- sum3x3Stencil :: Fractional a => Stencil Ix2 a a -- sum3x3Stencil = makeConvolutionStencil (3 :. 3) (1 :. 1) $ \ get -> @@ -50,7 +51,6 @@ i <- getDim sz r setDim zeroIndex r i - stencilSpec :: Spec stencilSpec = do describe "MapSingletonStencil" $ do @@ -78,8 +78,80 @@ (Default a, Unbox a, Manifest r Ix2 a) => Ix2 -> Ix2 -> Array r Ix2 a -> Array U Ix2 a stencilCorners ixC ix = computeAs U . mapStencil (Fill def) (makeStencil (3 :. 3) ixC $ \f -> f ix) +stencilConvolution :: Spec +stencilConvolution = do + let xs3 :: Array U Ix1 Int + xs3 = [1, 2, 3] + xs3f f = f (-1) 1 . f 0 2 . f 1 3 + xs4 :: Array U Ix1 Int + xs4 = [1, 2, 3, 4] + xs4f f = f (-2) 1 . f (-1) 2 . f 0 3 . f 1 4 + ys :: Array U Ix1 Int + ys = [1, 2, 3, 4, 5] + ysConvXs3 = [4, 10, 16, 22, 22] + ysConvXs4 = [10, 20, 30, 34, 31] + ysCorrXs3 = [8, 14, 20, 26, 14] + ysCorrXs4 = [11, 20, 30, 40, 26] + ysConvXs4' = [4, 10, 20, 30, 34] + ysCorrXs4' = [20, 30, 40, 26, 14] + xs4f' f = f (-1) 1 . f 0 2 . f 1 3 . f 2 4 + applyStencil s = computeAs U . mapStencil (Fill 0) s + describe "makeConvolutionStencilFromKernel" $ do + it "1x3" $ applyStencil (makeConvolutionStencilFromKernel xs3) ys `shouldBe` ysConvXs3 + it "1x4" $ applyStencil (makeConvolutionStencilFromKernel xs4) ys `shouldBe` ysConvXs4 + describe "makeCorrelationStencilFromKernel" $ do + it "1x3" $ applyStencil (makeCorrelationStencilFromKernel xs3) ys `shouldBe` ysCorrXs3 + it "1x4" $ applyStencil (makeCorrelationStencilFromKernel xs4) ys `shouldBe` ysCorrXs4 + describe "makeConvolutionStencil" $ do + it "1x3" $ applyStencil (makeConvolutionStencil (Sz1 3) 1 xs3f) ys `shouldBe` ysConvXs3 + it "1x4" $ applyStencil (makeConvolutionStencil (Sz1 4) 2 xs4f) ys `shouldBe` ysConvXs4 + it "1x4" $ applyStencil (makeConvolutionStencil (Sz1 4) 1 xs4f') ys `shouldBe` ysConvXs4' + describe "makeCorrelationStencil" $ do + it "1x3" $ applyStencil (makeCorrelationStencil (Sz1 3) 1 xs3f) ys `shouldBe` ysCorrXs3 + it "1x4" $ applyStencil (makeCorrelationStencil (Sz1 4) 2 xs4f) ys `shouldBe` ysCorrXs4 + it "1x4" $ applyStencil (makeCorrelationStencil (Sz1 4) 1 xs4f') ys `shouldBe` ysCorrXs4' + describe "makeConvolutionStencil == makeConvolutionStencilFromKernel" $ do + it "Sobel Horizontal" $ + property $ \(arr :: Array U Ix2 Int) -> + applyStencil (makeConvolutionStencil (3 :. 3) 1 sobelX) arr === + applyStencil (makeConvolutionStencilFromKernel sobelKernelX) arr + it "1x3" $ + property $ \(arr :: Array U Ix1 Int) -> + applyStencil (makeConvolutionStencil (Sz1 3) 1 xs3f) arr === + applyStencil (makeConvolutionStencilFromKernel xs3) arr + it "1x4" $ + property $ \(arr :: Array U Ix1 Int) -> + applyStencil (makeConvolutionStencil (Sz1 4) 2 xs4f) arr === + applyStencil (makeConvolutionStencilFromKernel xs4) arr + describe "makeCorrelationStencil == makeCorrelationStencilFromKernel" $ do + it "Sobel Horizontal" $ + property $ \(arr :: Array U Ix2 Int) -> + applyStencil (makeCorrelationStencil (3 :. 3) 1 sobelX) arr === + applyStencil (makeCorrelationStencilFromKernel sobelKernelX) arr + it "1x3" $ + property $ \(arr :: Array U Ix1 Int) -> + applyStencil (makeCorrelationStencil (Sz1 3) 1 xs3f) arr === + applyStencil (makeCorrelationStencilFromKernel xs3) arr + it "1x4" $ + property $ \(arr :: Array U Ix1 Int) -> + applyStencil (makeCorrelationStencil (Sz1 4) 2 xs4f) arr === + applyStencil (makeCorrelationStencilFromKernel xs4) arr + describe "makeConvolutionStencil == makeCorrelationStencil . rotate180" $ do + it "Sobel Horizontal" $ + property $ \(arr :: Array U Ix2 Int) -> + applyStencil (makeConvolutionStencilFromKernel sobelKernelX) arr === + applyStencil (makeCorrelationStencilFromKernel (rotate180 sobelKernelX)) arr + it "1x3" $ + property $ \(arr :: Array U Ix1 Int) -> + applyStencil (makeConvolutionStencilFromKernel xs3) arr === + applyStencil (makeCorrelationStencilFromKernel (rotate180 xs3)) arr + -- it "1x4" $ + -- property $ \(arr :: Array U Ix1 Int) -> + -- applyStencil (makeConvolutionStencilFromKernel xs4) arr === + -- applyStencil (makeCorrelationStencilFromKernel (rotate180 xs4)) arr + spec :: Spec -spec = +spec = do describe "Stencil" $ do stencilSpec let arr = [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: Array U Ix2 Int @@ -106,8 +178,24 @@ stride = Stride 2 it "map stencil with stride on small array" $ let strideArr = mapStencil (Fill 0) stencil arr - in computeWithStrideAs U stride strideArr `shouldBe` [[-4, 8],[2, 14]] + in computeWithStrideAs U stride strideArr `shouldBe` [[-4, 8], [2, 14]] it "map stencil with stride on larger array" $ let largeArr = makeArrayR U Seq (5 :. 5) (succ . toLinearIndex (5 :. 5)) strideArr = mapStencil (Fill 0) stencil largeArr - in computeWithStrideAs U stride strideArr `shouldBe` [[-6, 1, 14], [-13, 9, 43], [4, 21, 44]] + in computeWithStrideAs U stride strideArr `shouldBe` + [[-6, 1, 14], [-13, 9, 43], [4, 21, 44]] + stencilConvolution + +sobelX :: Num e => (Ix2 -> e -> e -> e) -> e -> e +sobelX f = f (-1 :. -1) (-1) . f (-1 :. 1) 1 . + f ( 0 :. -1) (-2) . f ( 0 :. 1) 2 . + f ( 1 :. -1) (-1) . f ( 1 :. 1) 1 + +sobelKernelX :: Array U Ix2 Int +sobelKernelX = [ [-1, 0, 1] + , [-2, 0, 2] + , [-1, 0, 1] ] + +rotate180 :: (Num ix, Index ix) => Array U ix Int -> Array U ix Int +rotate180 arr = computeAs U $ backpermute sz (\ix -> sz - 1 - ix) arr + where sz = size arr