diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,16 @@ +# 0.4.3 + +* Addition of `applyStencil` and `Padding` with helper functions `noPadding` and `samePadding`. +* Addition of `foldlStencil`, `foldrStencil` and monoidal `foldStencil`. +* Addition of common generic stencils: `sumStencil`, `productStencil`, `avgStencil`, + `maxStencil`, `minStencil` and `idStencil`. +* Addition of `mapStencilUnsafe` for the brave. +* Improve compile time error reporting for invalid dimensions. +* Fix incorrect loading of `DW` arrays of dimension higher than 3 +* Addition of `foldOuterSlice`, `ifoldOuterSlice`, `foldInnerSlice` and + `ifoldInnerSlice`. Fix for [#56](https://github.com/lehins/massiv/issues/56) + + # 0.4.2 * Fix loading empty `DS` stream arrays of unknown size. Fix for [#83](https://github.com/lehins/massiv/issues/83). diff --git a/massiv.cabal b/massiv.cabal --- a/massiv.cabal +++ b/massiv.cabal @@ -1,5 +1,5 @@ name: massiv -version: 0.4.2.0 +version: 0.4.3.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.hs b/src/Data/Massiv/Array.hs --- a/src/Data/Massiv/Array.hs +++ b/src/Data/Massiv/Array.hs @@ -222,5 +222,3 @@ All folding is done in a row-major order. -} - - diff --git a/src/Data/Massiv/Array/Delayed/Push.hs b/src/Data/Massiv/Array/Delayed/Push.hs --- a/src/Data/Massiv/Array/Delayed/Push.hs +++ b/src/Data/Massiv/Array/Delayed/Push.hs @@ -28,7 +28,6 @@ import Control.Monad import Data.Massiv.Core.Common -import Data.Massiv.Core.Index.Internal (Sz(SafeSz)) import Prelude hiding (map, zipWith) import Control.Scheduler as S (traverse_) import Data.Foldable as F diff --git a/src/Data/Massiv/Array/Delayed/Windowed.hs b/src/Data/Massiv/Array/Delayed/Windowed.hs --- a/src/Data/Massiv/Array/Delayed/Windowed.hs +++ b/src/Data/Massiv/Array/Delayed/Windowed.hs @@ -156,13 +156,14 @@ , dwWindow = Just $! Window - { windowStart = unSz (Sz (liftIndex2 min wStart (liftIndex (subtract 1) sz))) - , windowSize = Sz (liftIndex2 min wSize (liftIndex2 (-) sz wStart)) + { windowStart = wStart' + , windowSize = Sz (liftIndex2 min wSize (liftIndex2 (-) sz wStart')) , windowIndex = wIndex , windowUnrollIx2 = wUnrollIx2 } } where + wStart' = unSz (Sz (liftIndex2 min wStart (liftIndex (subtract 1) sz))) Sz sz = size arr Window { windowStart = wStart , windowSize = Sz wSize @@ -389,29 +390,31 @@ !(headWindowSz, tailWindowSz) = unconsSz windowSize !curWindowEnd = curWindowStart + unSz headWindowSz !pageElements = totalElem lowerSize - loadLower !i = - let !lowerWindow = - Window - { windowStart = lowerWindowStart - , windowSize = tailWindowSz - , windowIndex = windowIndex . consDim i - , windowUnrollIx2 = windowUnrollIx2 - } - !lowerArr = - DWArray - { dwArray = DArray Seq lowerSourceSize (indexBorder . consDim i) - , dwWindow = Just lowerWindow - } - in loadArrayWithStrideM - scheduler - (Stride lowerStrideIx) - lowerSize - lowerArr - (\k -> uWrite (k + pageElements * (i `div` s))) + mkLowerWindow i = + Window + { windowStart = lowerWindowStart + , windowSize = tailWindowSz + , windowIndex = windowIndex . consDim i + , windowUnrollIx2 = windowUnrollIx2 + } + mkLowerArray mw i = + DWArray + {dwArray = DArray Seq lowerSourceSize (indexBorder . consDim i), dwWindow = ($ i) <$> mw} + loadLower mw !i = + loadArrayWithStrideM + scheduler + (Stride lowerStrideIx) + lowerSize + (mkLowerArray mw i) + (\k -> uWrite (k + pageElements * (i `div` s))) {-# NOINLINE loadLower #-} - loopM_ 0 (< headDim windowStart) (+ s) loadLower - loopM_ (strideStart (Stride s) curWindowStart) (< curWindowEnd) (+ s) loadLower - loopM_ (strideStart (Stride s) curWindowEnd) (< unSz headSourceSize) (+ s) loadLower + loopM_ 0 (< headDim windowStart) (+ s) (loadLower Nothing) + loopM_ + (strideStart (Stride s) curWindowStart) + (< curWindowEnd) + (+ s) + (loadLower (Just mkLowerWindow)) + loopM_ (strideStart (Stride s) curWindowEnd) (< unSz headSourceSize) (+ s) (loadLower Nothing) {-# INLINE loadArrayWithIxN #-} @@ -430,24 +433,22 @@ !windowEnd = liftIndex2 (+) windowStart (unSz windowSize) !(t, windowStartL) = unconsDim windowStart !pageElements = totalElem szL - loadLower !i = - let !lowerWindow = - Window - { windowStart = windowStartL - , windowSize = snd $ unconsSz windowSize - , windowIndex = windowIndex . consDim i - , windowUnrollIx2 = windowUnrollIx2 - } - !lowerArr = - DWArray - { dwArray = DArray Seq szL (indexBorder . consDim i) - , dwWindow = Just lowerWindow - } - in with $ loadArrayM trivialScheduler_ lowerArr (\k -> uWrite (k + pageElements * i)) + mkLowerWindow i = + Window + { windowStart = windowStartL + , windowSize = snd $ unconsSz windowSize + , windowIndex = windowIndex . consDim i + , windowUnrollIx2 = windowUnrollIx2 + } + mkLowerArray mw i = + DWArray {dwArray = DArray Seq szL (indexBorder . consDim i), dwWindow = ($ i) <$> mw} + loadLower mw !i = + with $ + loadArrayM trivialScheduler_ (mkLowerArray mw i) (\k -> uWrite (k + pageElements * i)) {-# NOINLINE loadLower #-} - loopM_ 0 (< headDim windowStart) (+ 1) loadLower - loopM_ t (< headDim windowEnd) (+ 1) loadLower - loopM_ (headDim windowEnd) (< unSz si) (+ 1) loadLower + loopM_ 0 (< headDim windowStart) (+ 1) (loadLower Nothing) + loopM_ t (< headDim windowEnd) (+ 1) (loadLower (Just mkLowerWindow)) + loopM_ (headDim windowEnd) (< unSz si) (+ 1) (loadLower Nothing) {-# INLINE loadWithIxN #-} diff --git a/src/Data/Massiv/Array/Numeric/Integral.hs b/src/Data/Massiv/Array/Numeric/Integral.hs --- a/src/Data/Massiv/Array/Numeric/Integral.hs +++ b/src/Data/Massiv/Array/Numeric/Integral.hs @@ -134,11 +134,11 @@ -> Array r ix e -- ^ Array with values of @f(x,y,..)@ that will be used as source for integration. -> Array M ix e integralApprox stencil d sz n arr = - extract' zeroIndex sz $ toManifest $ loop 1 (<= coerce (dimensions sz)) (+ 1) arr applyStencil + extract' zeroIndex sz $ toManifest $ loop 1 (<= coerce (dimensions sz)) (+ 1) arr integrateAlong where !dx = d / fromIntegral n - applyStencil dim = integrateWith (stencil dx) (Dim dim) n - {-# INLINE applyStencil #-} + integrateAlong dim = integrateWith (stencil dx) (Dim dim) n + {-# INLINE integrateAlong #-} {-# INLINE integralApprox #-} diff --git a/src/Data/Massiv/Array/Ops/Fold.hs b/src/Data/Massiv/Array/Ops/Fold.hs --- a/src/Data/Massiv/Array/Ops/Fold.hs +++ b/src/Data/Massiv/Array/Ops/Fold.hs @@ -21,6 +21,10 @@ , foldMono , ifoldSemi , foldSemi + , foldOuterSlice + , ifoldOuterSlice + , foldInnerSlice + , ifoldInnerSlice , minimumM , minimum' , maximumM @@ -42,16 +46,19 @@ , foldlInner , ifoldrInner , foldrInner - -- *** Type safe + , foldInner + -- *** Type safe within , ifoldlWithin , foldlWithin , ifoldrWithin , foldrWithin - -- *** Partial + , foldWithin + -- *** Partial within , ifoldlWithin' , foldlWithin' , ifoldrWithin' , foldrWithin' + , foldWithin' -- ** Sequential folds @@ -91,6 +98,7 @@ ) where import Data.Massiv.Array.Delayed.Pull +import Data.Massiv.Array.Ops.Construct import Data.Massiv.Array.Ops.Fold.Internal import Data.Massiv.Core import Data.Massiv.Core.Common @@ -281,6 +289,106 @@ foldrInner = foldrWithin' 1 {-# INLINE foldrInner #-} +-- | Monoidal fold over the inner most dimension. +-- +-- @since 0.4.3 +foldInner :: (Monoid e, Index (Lower ix), Source r ix e) => Array r ix e -> Array D (Lower ix) e +foldInner = foldlInner mappend mempty +{-# INLINE foldInner #-} + +-- | Monoidal fold over some internal dimension. +-- +-- @since 0.4.3 +foldWithin :: + (Source r ix a, Monoid a, Index (Lower ix), IsIndexDimension ix n) + => Dimension n + -> Array r ix a + -> Array D (Lower ix) a +foldWithin dim = foldlWithin dim mappend mempty +{-# INLINE foldWithin #-} + +-- | Monoidal fold over some internal dimension. This is a pratial function and will +-- result in `IndexDimensionException` if supplied dimension is invalid. +-- +-- @since 0.4.3 +foldWithin' :: + (Source r ix a, Monoid a, Index (Lower ix)) + => Dim + -> Array r ix a + -> Array D (Lower ix) a +foldWithin' dim = foldlWithin' dim mappend mempty +{-# INLINE foldWithin' #-} + + +-- | Reduce each outer slice into a monoid and mappend results together +-- +-- ==== __Example__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> import Data.Monoid (Product(..)) +-- >>> arr = computeAs P $ iterateN (Sz2 2 3) (+1) (10 :: Int) +-- >>> arr +-- Array P Seq (Sz (2 :. 3)) +-- [ [ 11, 12, 13 ] +-- , [ 14, 15, 16 ] +-- ] +-- >>> getProduct $ foldOuterSlice (\row -> Product (A.sum row)) arr +-- 1620 +-- >>> (11 + 12 + 13) * (14 + 15 + 16) :: Int +-- 1620 +-- +-- @since 0.4.3 +foldOuterSlice :: (OuterSlice r ix e, Monoid m) => (Elt r ix e -> m) -> Array r ix e -> m +foldOuterSlice f = ifoldOuterSlice (const f) +{-# INLINE foldOuterSlice #-} + + +-- | Reduce each outer slice into a monoid with an index aware function and mappend results +-- together +-- +-- @since 0.4.3 +ifoldOuterSlice :: (OuterSlice r ix e, Monoid m) => (Ix1 -> Elt r ix e -> m) -> Array r ix e -> m +ifoldOuterSlice f arr = foldMono g $ range (getComp arr) 0 (headDim (unSz (size arr))) + where + g i = f i (unsafeOuterSlice arr i) + {-# INLINE g #-} +{-# INLINE ifoldOuterSlice #-} + + +-- | Reduce each inner slice into a monoid and mappend results together +-- +-- ==== __Example__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> import Data.Monoid (Product(..)) +-- >>> arr = computeAs P $ iterateN (Sz2 2 3) (+1) (10 :: Int) +-- >>> arr +-- Array P Seq (Sz (2 :. 3)) +-- [ [ 11, 12, 13 ] +-- , [ 14, 15, 16 ] +-- ] +-- >>> getProduct $ foldInnerSlice (\column -> Product (A.sum column)) arr +-- 19575 +-- >>> (11 + 14) * (12 + 15) * (13 + 16) :: Int +-- 19575 +-- +-- @since 0.4.3 +foldInnerSlice :: (InnerSlice r ix e, Monoid m) => (Elt r ix e -> m) -> Array r ix e -> m +foldInnerSlice f = ifoldInnerSlice (const f) +{-# INLINE foldInnerSlice #-} + + +-- | Reduce each inner slice into a monoid with an index aware function and mappend +-- results together +-- +-- @since 0.4.3 +ifoldInnerSlice :: (InnerSlice r ix e, Monoid m) => (Ix1 -> Elt r ix e -> m) -> Array r ix e -> m +ifoldInnerSlice f arr = foldMono g $ range (getComp arr) 0 (unSz k) + where + szs@(_, !k) = unsnocSz (size arr) + g i = f i (unsafeInnerSlice arr szs i) + {-# INLINE g #-} +{-# INLINE ifoldInnerSlice #-} -- | /O(n)/ - Compute maximum of all elements. -- diff --git a/src/Data/Massiv/Array/Ops/Transform.hs b/src/Data/Massiv/Array/Ops/Transform.hs --- a/src/Data/Massiv/Array/Ops/Transform.hs +++ b/src/Data/Massiv/Array/Ops/Transform.hs @@ -39,8 +39,12 @@ -- ** Append/Split , cons , unconsM + -- , headM + -- , head' , snoc , unsnocM + -- , lastM + -- , last' , appendM , append' , concatM @@ -406,6 +410,40 @@ uWrite startAt e >> dlLoad arr scheduler (startAt + 1) uWrite } {-# INLINE cons #-} + + +-- -- | /O(1)/ - Take the first element off the vector from the left side. +-- -- +-- -- @since 0.4.3 +-- headM :: (MonadThrow m, Source r Ix1 e) => Array r Ix1 e -> m e +-- headM = fmap fst . unconsM +-- {-# INLINE headM #-} + +-- -- | /O(1)/ - Take the first element off the vector from the left side. Throws +-- -- `SizeEmptyException` +-- -- +-- -- @since 0.4.3 +-- head' :: Source r Ix1 e => Array r Ix1 e -> e +-- head' = either throw id . headM +-- {-# INLINE head' #-} + + +-- -- | /O(1)/ - Take the last element off the vector from the right side. +-- -- +-- -- @since 0.4.3 +-- lastM :: (MonadThrow m, Source r Ix1 e) => Array r Ix1 e -> m e +-- lastM = fmap snd . unsnocM +-- {-# INLINE lastM #-} + +-- -- | /O(1)/ - Take the last element off the vector from the right side. Throws +-- -- `SizeEmptyException` +-- -- +-- -- @since 0.4.3 +-- last' :: Source r Ix1 e => Array r Ix1 e -> e +-- last' = either throw id . lastM +-- {-# INLINE last' #-} + + -- | /O(1)/ - Take one element off the vector from the left side. -- 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 @@ -16,8 +16,26 @@ , Value , makeStencil , makeStencilDef + , getStencilSize + , getStencilCenter + -- ** Padding + , Padding(..) + , noPadding + , samePadding + -- ** Application , mapStencil - -- ** Profunctor + , applyStencil + -- ** Common stencils + , idStencil + , sumStencil + , productStencil + , avgStencil + , maxStencil + , minStencil + , foldlStencil + , foldrStencil + , foldStencil + -- ** Profunctor , dimapStencil , lmapStencil , rmapStencil @@ -27,18 +45,32 @@ , Default(def) ) where +import Data.Coerce import Data.Default.Class (Default(def)) import Data.Massiv.Array.Delayed.Windowed import Data.Massiv.Array.Manifest import Data.Massiv.Array.Stencil.Convolution import Data.Massiv.Array.Stencil.Internal +import Data.Massiv.Array.Stencil.Unsafe import Data.Massiv.Core.Common +import Data.Semigroup import GHC.Exts (inline) +-- | Get the size of the stencil +-- +-- @since 0.4.3 +getStencilSize :: Stencil ix e a -> Sz ix +getStencilSize = stencilSize --- | Map a constructed stencil over an array. Resulting array must be `compute`d in order to be --- useful. +-- | Get the index of the stencil's center -- +-- @since 0.4.3 +getStencilCenter :: Stencil ix e a -> ix +getStencilCenter = stencilCenter + +-- | Map a constructed stencil over an array. Resulting array must be +-- `Data.Massiv.Array.compute`d in order to be useful. +-- -- @since 0.1.0 mapStencil :: (Source r ix e, Manifest r ix e) @@ -46,19 +78,137 @@ -> Stencil ix e a -- ^ Stencil to map over the array -> Array r ix e -- ^ Source array -> Array DW ix a -mapStencil b (Stencil sSz sCenter stencilF) !arr = insertWindow warr window +mapStencil b stencil = applyStencil (samePadding stencil b) stencil +{-# INLINE mapStencil #-} + + +-- | Padding of the source array before stencil application. +-- +-- ==== __Examples__ +-- +-- In order to see the affect of padding we can simply apply an identity stencil to an +-- array: +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ resize' (Sz2 2 3) (Ix1 1 ... 6) +-- >>> applyStencil noPadding idStencil a +-- Array DW Seq (Sz (2 :. 3)) +-- [ [ 1, 2, 3 ] +-- , [ 4, 5, 6 ] +-- ] +-- >>> applyStencil (Padding (Sz2 1 2) (Sz2 3 4) (Fill 0)) idStencil a +-- Array DW Seq (Sz (6 :. 9)) +-- [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ] +-- , [ 0, 0, 1, 2, 3, 0, 0, 0, 0 ] +-- , [ 0, 0, 4, 5, 6, 0, 0, 0, 0 ] +-- , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ] +-- , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ] +-- , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ] +-- ] +-- +-- It is also a nice technique to see the border resolution strategy in action: +-- +-- >>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Wrap) idStencil a +-- Array DW Seq (Sz (6 :. 9)) +-- [ [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ] +-- , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ] +-- , [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ] +-- , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ] +-- , [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ] +-- , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ] +-- ] +-- >>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Edge) idStencil a +-- Array DW Seq (Sz (6 :. 9)) +-- [ [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ] +-- , [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ] +-- , [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ] +-- , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ] +-- , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ] +-- , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ] +-- ] +-- >>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Reflect) idStencil a +-- Array DW Seq (Sz (6 :. 9)) +-- [ [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ] +-- , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ] +-- , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ] +-- , [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ] +-- , [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ] +-- , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ] +-- ] +-- >>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Continue) idStencil a +-- Array DW Seq (Sz (6 :. 9)) +-- [ [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ] +-- , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ] +-- , [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ] +-- , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ] +-- , [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ] +-- , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ] +-- ] +-- +-- @since 0.4.3 +data Padding ix e = Padding + { paddingFromOrigin :: !(Sz ix) + , paddingFromBottom :: !(Sz ix) + , paddingWithElement :: !(Border e) + -- ^ Element to do padding with + } deriving (Eq, Show) + +-- | Also known as "valid" padding. When stencil is applied to an array, that array will +-- shrink, unless the stencil is of size 1. +-- +-- @since 0.4.3 +noPadding :: Index ix => Padding ix e +noPadding = Padding zeroSz zeroSz Edge + +-- | Padding that matches the size of the stencil, which is known as "same" padding, +-- because when a stencil is applied to an array with such matching padding, the resulting +-- array will be of the same size as the source array. This is exactly the behavior of +-- `mapStencil` +-- +-- @since 0.4.3 +samePadding :: Index ix => Stencil ix e a -> Border e -> Padding ix e +samePadding (Stencil (Sz sSz) sCenter _) border = + Padding + { paddingFromOrigin = Sz sCenter + , paddingFromBottom = Sz (liftIndex2 (-) sSz (liftIndex (+1) sCenter)) + , paddingWithElement = border + } + +-- | Apply a constructed stencil over an array. Resulting array must be +-- `Data.Massiv.Array.compute`d in order to be useful. Unlike `mapStencil`, the size of +-- the resulting array will not necesserally be the same as the source array, which will +-- depend on the padding. +-- +-- @since 0.4.3 +applyStencil :: + (Source r ix e, Manifest r ix e) + => Padding ix e + -- ^ Padding to be applied to the source array. This will dictate the resulting size of + -- the array. No padding will cause it to shrink by the size of the stencil + -> Stencil ix e a -- ^ Stencil to apply to the array + -> Array r ix e -- ^ Source array + -> Array DW ix a +applyStencil (Padding (Sz po) (Sz pb) border) (Stencil sSz sCenter stencilF) !arr = + insertWindow warr window where - !warr = DArray (getComp arr) sz (unValue . stencilF (Value . borderIndex b arr)) + !offset = liftIndex2 (-) sCenter po + !warr = + DArray + (getComp arr) + sz + (unValue . stencilF (Value . borderIndex border arr) . liftIndex2 (+) offset) + -- Size by which the resulting array will shrink (not accounting for padding) + !shrinkSz = Sz (liftIndex (subtract 1) (unSz sSz)) + !sz = liftSz2 (-) (SafeSz (liftIndex2 (+) po (liftIndex2 (+) pb (unSz (size arr))))) shrinkSz + !wsz = liftSz2 (-) (size arr) shrinkSz !window = Window - { windowStart = sCenter - , windowSize = windowSz - , windowIndex = unValue . stencilF (Value . unsafeIndex arr) + { windowStart = po + , windowSize = wsz + , windowIndex = unValue . stencilF (Value . unsafeIndex arr) . liftIndex2 (+) offset , windowUnrollIx2 = unSz . fst <$> pullOutSzM sSz 2 } - !sz = size arr - !windowSz = Sz (liftIndex2 (-) (unSz sz) (liftIndex (subtract 1) (unSz sSz))) -{-# INLINE mapStencil #-} +{-# INLINE applyStencil #-} -- | Construct a stencil from a function, which describes how to calculate the @@ -71,9 +221,10 @@ -- -- Below is an example of creating a `Stencil`, which, when mapped over a -- 2-dimensional array, will compute an average of all elements in a 3x3 square --- for each element in that array. /Note:/ Make sure to add @INLINE@ pragma, --- otherwise performance will be terrible. +-- for each element in that array. -- +-- /Note/ - Make sure to add an @INLINE@ pragma, otherwise performance will be terrible. +-- -- > average3x3Stencil :: (Default a, Fractional a) => Stencil Ix2 a a -- > average3x3Stencil = makeStencil (Sz (3 :. 3)) (1 :. 1) $ \ get -> -- > ( get (-1 :. -1) + get (-1 :. 0) + get (-1 :. 1) + @@ -102,7 +253,7 @@ -- @since 0.2.3 makeStencilDef :: Index ix - => e + => e -- ^ Default element that will be used for stencil validation only. -> Sz ix -- ^ Size of the stencil -> ix -- ^ Center of the stencil -> ((ix -> Value e) -> Value a) @@ -116,3 +267,218 @@ {-# INLINE stencil #-} {-# INLINE makeStencilDef #-} +-- | Identity stencil that does not change the elements of the source array. +-- +-- @since 0.4.3 +idStencil :: Index ix => Stencil ix e e +idStencil = makeUnsafeStencil oneSz zeroIndex $ \ _ get -> get zeroIndex +{-# INLINE idStencil #-} + +-- | Stencil that does a left fold in a row-major order. Regardless of the supplied size +-- resulting stencil will be centered at zero, although by using `Padding` it is possible +-- to overcome this limitation. +-- +-- ==== __Examples__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Int) +-- >>> a +-- Array P Seq (Sz (3 :. 4)) +-- [ [ 11, 12, 13, 14 ] +-- , [ 15, 16, 17, 18 ] +-- , [ 19, 20, 21, 22 ] +-- ] +-- >>> applyStencil noPadding (foldlStencil (flip (:)) [] (Sz2 3 2)) a +-- Array DW Seq (Sz (1 :. 3)) +-- [ [ [20,19,16,15,12,11], [21,20,17,16,13,12], [22,21,18,17,14,13] ] +-- ] +-- >>> applyStencil (Padding (Sz2 1 0) 0 (Fill 10)) (foldlStencil (flip (:)) [] (Sz2 3 2)) a +-- Array DW Seq (Sz (2 :. 3)) +-- [ [ [16,15,12,11,10,10], [17,16,13,12,10,10], [18,17,14,13,10,10] ] +-- , [ [20,19,16,15,12,11], [21,20,17,16,13,12], [22,21,18,17,14,13] ] +-- ] +-- +-- @since 0.4.3 +foldlStencil :: Index ix => (a -> e -> a) -> a -> Sz ix -> Stencil ix e a +foldlStencil f acc0 sz = + makeUnsafeStencil sz zeroIndex $ \_ get -> + iter zeroIndex (unSz sz) oneIndex (<) acc0 $ \ix -> (`f` get ix) +{-# INLINE foldlStencil #-} + +-- | Stencil that does a right fold in a row-major order. Regardless of the supplied size +-- resulting stencil will be centered at zero, although by using `Padding` it is possible +-- to overcome this limitation. +-- +-- ==== __Examples__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Int) +-- >>> a +-- Array P Seq (Sz (3 :. 4)) +-- [ [ 11, 12, 13, 14 ] +-- , [ 15, 16, 17, 18 ] +-- , [ 19, 20, 21, 22 ] +-- ] +-- >>> applyStencil noPadding (foldrStencil (:) [] (Sz2 2 3)) a +-- Array DW Seq (Sz (2 :. 2)) +-- [ [ [11,12,13,15,16,17], [12,13,14,16,17,18] ] +-- , [ [15,16,17,19,20,21], [16,17,18,20,21,22] ] +-- ] +-- +-- @since 0.4.3 +foldrStencil :: Index ix => (e -> a -> a) -> a -> Sz ix -> Stencil ix e a +foldrStencil f acc0 sz = + let ixStart = liftIndex2 (-) (unSz sz) oneIndex + in makeUnsafeStencil sz zeroIndex $ \_ get -> + iter ixStart zeroIndex (pureIndex (-1)) (>=) acc0 $ \ix -> f (get ix) +{-# INLINE foldrStencil #-} + + +foldStencil :: (Monoid e, Index ix) => Sz ix -> Stencil ix e e +foldStencil = foldlStencil mappend mempty +{-# INLINE foldStencil #-} + +-- | Create a stencil centered at 0 that will extract the maximum value in the region of +-- supplied size. +-- +-- ==== __Example__ +-- +-- Here is a sample implementation of min pooling. +-- +-- >>> import Data.Massiv.Array as A +-- >>> a <- computeAs P <$> resizeM (Sz2 9 9) (Ix1 10 ..: 91) +-- >>> a +-- Array P Seq (Sz (9 :. 9)) +-- [ [ 10, 11, 12, 13, 14, 15, 16, 17, 18 ] +-- , [ 19, 20, 21, 22, 23, 24, 25, 26, 27 ] +-- , [ 28, 29, 30, 31, 32, 33, 34, 35, 36 ] +-- , [ 37, 38, 39, 40, 41, 42, 43, 44, 45 ] +-- , [ 46, 47, 48, 49, 50, 51, 52, 53, 54 ] +-- , [ 55, 56, 57, 58, 59, 60, 61, 62, 63 ] +-- , [ 64, 65, 66, 67, 68, 69, 70, 71, 72 ] +-- , [ 73, 74, 75, 76, 77, 78, 79, 80, 81 ] +-- , [ 82, 83, 84, 85, 86, 87, 88, 89, 90 ] +-- ] +-- >>> computeWithStrideAs P (Stride 3) $ mapStencil Edge (maxStencil (Sz 3)) a +-- Array P Seq (Sz (3 :. 3)) +-- [ [ 30, 33, 36 ] +-- , [ 57, 60, 63 ] +-- , [ 84, 87, 90 ] +-- ] +-- +-- @since 0.4.3 +maxStencil :: (Bounded e, Ord e, Index ix) => Sz ix -> Stencil ix e e +maxStencil = dimapStencil coerce getMax . foldStencil +{-# INLINE maxStencil #-} + + +-- | Create a stencil centered at 0 that will extract the maximum value in the region of +-- supplied size. +-- +-- ==== __Example__ +-- +-- Here is a sample implementation of min pooling. +-- +-- >>> import Data.Massiv.Array as A +-- >>> a <- computeAs P <$> resizeM (Sz2 9 9) (Ix1 10 ..: 91) +-- >>> a +-- Array P Seq (Sz (9 :. 9)) +-- [ [ 10, 11, 12, 13, 14, 15, 16, 17, 18 ] +-- , [ 19, 20, 21, 22, 23, 24, 25, 26, 27 ] +-- , [ 28, 29, 30, 31, 32, 33, 34, 35, 36 ] +-- , [ 37, 38, 39, 40, 41, 42, 43, 44, 45 ] +-- , [ 46, 47, 48, 49, 50, 51, 52, 53, 54 ] +-- , [ 55, 56, 57, 58, 59, 60, 61, 62, 63 ] +-- , [ 64, 65, 66, 67, 68, 69, 70, 71, 72 ] +-- , [ 73, 74, 75, 76, 77, 78, 79, 80, 81 ] +-- , [ 82, 83, 84, 85, 86, 87, 88, 89, 90 ] +-- ] +-- >>> computeWithStrideAs P (Stride 3) $ mapStencil Edge (minStencil (Sz 3)) a +-- Array P Seq (Sz (3 :. 3)) +-- [ [ 10, 13, 16 ] +-- , [ 37, 40, 43 ] +-- , [ 64, 67, 70 ] +-- ] +-- +-- @since 0.4.3 +minStencil :: (Bounded e, Ord e, Index ix) => Sz ix -> Stencil ix e e +minStencil = dimapStencil coerce getMin . foldStencil +{-# INLINE minStencil #-} + +-- | Sum all elements in the stencil region +-- +-- ==== __Examples__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ iterateN (Sz2 2 5) (* 2) (1 :: Int) +-- >>> a +-- Array P Seq (Sz (2 :. 5)) +-- [ [ 2, 4, 8, 16, 32 ] +-- , [ 64, 128, 256, 512, 1024 ] +-- ] +-- >>> applyStencil noPadding (sumStencil (Sz2 1 2)) a +-- Array DW Seq (Sz (2 :. 4)) +-- [ [ 6, 12, 24, 48 ] +-- , [ 192, 384, 768, 1536 ] +-- ] +-- >>> [2 + 4, 4 + 8, 8 + 16, 16 + 32] :: [Int] +-- [6,12,24,48] +-- +-- @since 0.4.3 +sumStencil :: (Num e, Index ix) => Sz ix -> Stencil ix e e +sumStencil = dimapStencil coerce getSum . foldStencil +{-# INLINE sumStencil #-} + +-- | Multiply all elements in the stencil region +-- +-- ==== __Examples__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ iterateN (Sz2 2 2) (+1) (0 :: Int) +-- >>> a +-- Array P Seq (Sz (2 :. 2)) +-- [ [ 1, 2 ] +-- , [ 3, 4 ] +-- ] +-- >>> applyStencil (Padding 0 2 (Fill 0)) (productStencil 2) a +-- Array DW Seq (Sz (3 :. 3)) +-- [ [ 24, 0, 0 ] +-- , [ 0, 0, 0 ] +-- , [ 0, 0, 0 ] +-- ] +-- >>> applyStencil (Padding 0 2 Reflect) (productStencil 2) a +-- Array DW Seq (Sz (3 :. 3)) +-- [ [ 24, 64, 24 ] +-- , [ 144, 256, 144 ] +-- , [ 24, 64, 24 ] +-- ] +-- +-- @since 0.4.3 +productStencil :: (Num e, Index ix) => Sz ix -> Stencil ix e e +productStencil = dimapStencil coerce getProduct . foldStencil +{-# INLINE productStencil #-} + +-- | Find the average value of all elements in the stencil region +-- +-- ==== __Example__ +-- +-- >>> import Data.Massiv.Array as A +-- >>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Double) +-- >>> a +-- Array P Seq (Sz (3 :. 4)) +-- [ [ 11.0, 12.0, 13.0, 14.0 ] +-- , [ 15.0, 16.0, 17.0, 18.0 ] +-- , [ 19.0, 20.0, 21.0, 22.0 ] +-- ] +-- >>> applyStencil noPadding (avgStencil (Sz2 2 3)) a +-- Array DW Seq (Sz (2 :. 2)) +-- [ [ 14.0, 15.0 ] +-- , [ 18.0, 19.0 ] +-- ] +-- >>> Prelude.sum [11.0, 12.0, 13.0, 15.0, 16.0, 17.0] / 6 :: Double +-- 14.0 +-- +-- @since 0.4.3 +avgStencil :: (Fractional e, Index ix) => Sz ix -> Stencil ix e e +avgStencil sz = sumStencil sz / fromIntegral (totalElem sz) +{-# INLINE avgStencil #-} diff --git a/src/Data/Massiv/Array/Stencil/Internal.hs b/src/Data/Massiv/Array/Stencil/Internal.hs --- a/src/Data/Massiv/Array/Stencil/Internal.hs +++ b/src/Data/Massiv/Array/Stencil/Internal.hs @@ -269,7 +269,9 @@ validateStencil :: Index ix => e -> Stencil ix e a -> Stencil ix e a -validateStencil d s@(Stencil sSz sCenter stencil) = - let valArr = DArray Seq sSz (const d) - in stencil (Value . safeStencilIndex valArr) sCenter `seq` s +validateStencil d s@(Stencil sSz sCenter stencil) + | isSafeIndex sSz sCenter = + let valArr = DArray Seq sSz (const d) + in stencil (Value . safeStencilIndex valArr) sCenter `seq` s + | otherwise = throw $ IndexOutOfBoundsException sSz sCenter {-# INLINE validateStencil #-} 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 @@ -13,6 +13,7 @@ module Data.Massiv.Array.Stencil.Unsafe ( -- * Stencil makeUnsafeStencil + , mapStencilUnsafe , forStencilUnsafe ) where @@ -23,9 +24,9 @@ import GHC.Exts (inline) --- | This is an unsafe version of the stencil computation. There are no bounds checking further from --- the border, so if you do make sure you are not going outside the size of the stencil, you will be --- safe, but this is not enforced. +-- | Just as `mapStencilUnsafe` this is an unsafe version of the stencil +-- mapping. Arguments are in slightly different order and the indexing function returns +-- `Nothing` for elements outside the border. -- -- @since 0.1.7 forStencilUnsafe :: @@ -56,6 +57,35 @@ {-# INLINE stencil #-} {-# INLINE forStencilUnsafe #-} + +-- | This is an unsafe version of `Data.Massiv.Array.Stencil.mapStencil`, that does no +-- stencil validation. There is no performance difference between the two, but the unsafe +-- version has an advantage of not requiring to deal with `Value` wrapper. +-- +-- @since 0.4.3 +mapStencilUnsafe :: + Manifest r ix e + => Border e + -> Sz ix + -> ix + -> ((ix -> e) -> a) + -> Array r ix e + -> Array DW ix a +mapStencilUnsafe b sSz sCenter stencilF !arr = insertWindow warr window + where + !warr = DArray (getComp arr) sz (stencil (borderIndex b arr)) + !window = + Window + { windowStart = sCenter + , windowSize = windowSz + , windowIndex = stencil (unsafeIndex arr) + , windowUnrollIx2 = unSz . fst <$> pullOutSzM sSz 2 + } + !sz = size arr + !windowSz = Sz (liftIndex2 (-) (unSz sz) (liftIndex (subtract 1) (unSz sSz))) + stencil getVal !ix = inline stencilF $ \ !ixD -> getVal (liftIndex2 (+) ix ixD) + {-# INLINE stencil #-} +{-# INLINE mapStencilUnsafe #-} -- | Similar to `Data.Massiv.Array.Stencil.makeStencil`, but there are no guarantees that the diff --git a/src/Data/Massiv/Core.hs b/src/Data/Massiv/Core.hs --- a/src/Data/Massiv/Core.hs +++ b/src/Data/Massiv/Core.hs @@ -47,7 +47,7 @@ ) where import Control.Exception (Exception(..), SomeException) -import Control.Scheduler (WorkerStates, initWorkerStates) +import Control.Scheduler (initWorkerStates) import Data.Massiv.Core.Common import Data.Massiv.Core.Index import Data.Massiv.Core.List diff --git a/src/Data/Massiv/Core/Index.hs b/src/Data/Massiv/Core/Index.hs --- a/src/Data/Massiv/Core/Index.hs +++ b/src/Data/Massiv/Core/Index.hs @@ -32,6 +32,7 @@ , zeroSz , oneSz , liftSz + , liftSz2 , consSz , unconsSz , snocSz @@ -43,6 +44,8 @@ , Dim(..) , Dimension(Dim1, Dim2, Dim3, Dim4, Dim5, DimN) , IsIndexDimension + , IsDimValid + , ReportInvalidDim -- ** Stride , Stride(Stride) , unStride @@ -480,7 +483,7 @@ => ix -- ^ Start index -> ix -- ^ End index -> ix -- ^ Increment - -> (Int -> Int -> Bool) -- ^ Continuation confition + -> (Int -> Int -> Bool) -- ^ Continuation condition -> a -- ^ Accumulator -> (ix -> a -> a) -- ^ Iterating function -> a diff --git a/src/Data/Massiv/Core/Index/Internal.hs b/src/Data/Massiv/Core/Index/Internal.hs --- a/src/Data/Massiv/Core/Index/Internal.hs +++ b/src/Data/Massiv/Core/Index/Internal.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE CPP #-} {-# LANGUAGE ConstraintKinds #-} @@ -12,7 +13,7 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeSynonymInstances #-} - +{-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-} #if __GLASGOW_HASKELL__ < 820 {-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} #endif @@ -33,6 +34,7 @@ , zeroSz , oneSz , liftSz + , liftSz2 , consSz , unconsSz , snocSz @@ -48,6 +50,8 @@ , pattern Dim4 , pattern Dim5 , IsIndexDimension + , IsDimValid + , ReportInvalidDim , Lower , Index(..) , Ix0(..) @@ -60,7 +64,7 @@ ) where import Control.DeepSeq -import Control.Exception (Exception(..)) +import Control.Exception (Exception(..), throw) import Control.Monad.Catch (MonadThrow(..)) import Data.Coerce import Data.Massiv.Core.Iterator @@ -182,7 +186,20 @@ liftSz f (SafeSz ix) = Sz (liftIndex f ix) {-# INLINE liftSz #-} +-- | Same as `liftIndex2`, but for `Sz` +-- +-- ==== __Example__ +-- +-- >>> import Data.Massiv.Core.Index +-- >>> liftSz2 (-) (Sz2 2 3) (Sz2 3 1) +-- Sz (0 :. 2) +-- +-- @since 0.4.3 +liftSz2 :: Index ix => (Int -> Int -> Int) -> Sz ix -> Sz ix -> Sz ix +liftSz2 f sz1 sz2 = Sz (liftIndex2 f (coerce sz1) (coerce sz2)) +{-# INLINE liftSz2 #-} + -- | Same as `consDim`, but for `Sz` -- -- ==== __Example__ @@ -329,7 +346,7 @@ -- safely used with it. -- -- @since 0.2.4 -type IsIndexDimension ix n = (1 <= n, n <= Dimensions ix, Index ix, KnownNat n) +type IsIndexDimension ix n = (IsDimValid ix n ~ 'True, Index ix, KnownNat n) -- | This type family will always point to a type for a dimension that is one lower than the type @@ -338,6 +355,18 @@ -- @since 0.1.0 type family Lower ix :: * + +type family ReportInvalidDim (dims :: Nat) (n :: Nat) isNotZero isLess :: Bool where + ReportInvalidDim dims n True True = True + ReportInvalidDim dims n True False = + TypeError (Text "Dimension " :<>: ShowType n :<>: Text " is higher than " :<>: + Text "the maximum expected " :<>: ShowType dims) + ReportInvalidDim dims n False isLess = + TypeError (Text "Zero dimensional indices are not supported") + +type family IsDimValid ix n :: Bool where + IsDimValid ix n = ReportInvalidDim (Dimensions ix) n (1 <=? n) (n <=? Dimensions ix) + -- | This is bread and butter of multi-dimensional array indexing. It is unlikely that any of the -- functions in this class will be useful to a regular user, unless general algorithms are being -- implemented that do span multiple dimensions. @@ -622,6 +651,7 @@ fromLinearIndexAcc n k = k `quotRem` n {-# INLINE [1] fromLinearIndexAcc #-} repairIndex k@(SafeSz ksz) !i rBelow rOver + | ksz <= 0 = throw $ IndexZeroException ksz | i < 0 = rBelow k i | i >= ksz = rOver k i | otherwise = i