diff --git a/Data/Array/Repa/Algorithms/ColorRamp.hs b/Data/Array/Repa/Algorithms/ColorRamp.hs --- a/Data/Array/Repa/Algorithms/ColorRamp.hs +++ b/Data/Array/Repa/Algorithms/ColorRamp.hs @@ -1,49 +1,49 @@ {-# LANGUAGE RankNTypes #-} -- | Hyprometric color ramps, for making pretty images from scalar data. -module Data.Array.Repa.Algorithms.ColorRamp - (rampColorHotToCold) +module Data.Array.Repa.Algorithms.ColorRamp + (rampColorHotToCold) where -- | Standard Hot to Cold hypsometric color ramp. --- Color sequence is red, yellow, green, cyan, blue. +-- Color sequence is red, yellow, green, cyan, blue. rampColorHotToCold - :: forall a - . (Ord a, Floating a) - => a -- ^ Minimum value of range. - -> a -- ^ Maximum value of range. - -> a -- ^ Data value. - -> (a, a, a) - + :: forall a + . (Ord a, Floating a) + => a -- ^ Minimum value of range. + -> a -- ^ Maximum value of range. + -> a -- ^ Data value. + -> (a, a, a) + {-# INLINE rampColorHotToCold #-} rampColorHotToCold vmin vmax vNotNorm - = let - v | vNotNorm < vmin = vmin - | vNotNorm > vmax = vmax - | otherwise = vNotNorm - - dv = vmax - vmin + = let + v | vNotNorm < vmin = vmin + | vNotNorm > vmax = vmax + | otherwise = vNotNorm + + dv = vmax - vmin - result | v < vmin + 0.25 * dv - = ( 0 - , 4 * (v - vmin) / dv - , 1.0) - - | v < vmin + 0.5 * dv - = ( 0 - , 1.0 - , 1 + 4 * (vmin + 0.25 * dv - v) / dv) - - | v < vmin + 0.75 * dv - = ( 4 * (v - vmin - 0.5 * dv) / dv - , 1.0 - , 0.0) - - | otherwise - = ( 1.0 - , 1 + 4 * (vmin + 0.75 * dv - v) / dv - , 0) - - in result + result | v < vmin + 0.25 * dv + = ( 0 + , 4 * (v - vmin) / dv + , 1.0) + + | v < vmin + 0.5 * dv + = ( 0 + , 1.0 + , 1 + 4 * (vmin + 0.25 * dv - v) / dv) + + | v < vmin + 0.75 * dv + = ( 4 * (v - vmin - 0.5 * dv) / dv + , 1.0 + , 0.0) + + | otherwise + = ( 1.0 + , 1 + 4 * (vmin + 0.75 * dv - v) / dv + , 0) + + in result diff --git a/Data/Array/Repa/Algorithms/Complex.hs b/Data/Array/Repa/Algorithms/Complex.hs --- a/Data/Array/Repa/Algorithms/Complex.hs +++ b/Data/Array/Repa/Algorithms/Complex.hs @@ -2,51 +2,51 @@ {-# OPTIONS -fno-warn-orphans #-} -- | Strict complex doubles. module Data.Array.Repa.Algorithms.Complex - ( Complex - , mag - , arg) + ( Complex + , mag + , arg) where -- | Complex doubles. type Complex - = (Double, Double) + = (Double, Double) instance Num Complex where {-# INLINE abs #-} - abs x = (mag x, 0) + abs x = (mag x, 0) {-# INLINE signum #-} - signum (re, _) = (signum re, 0) + signum (re, _) = (signum re, 0) {-# INLINE fromInteger #-} - fromInteger n = (fromInteger n, 0.0) + fromInteger n = (fromInteger n, 0.0) {-# INLINE (+) #-} - (r, i) + (r', i') = (r+r', i+i') + (r, i) + (r', i') = (r+r', i+i') {-# INLINE (-) #-} - (r, i) - (r', i') = (r-r', i-i') + (r, i) - (r', i') = (r-r', i-i') {-# INLINE (*) #-} - (r, i) * (r', i') = (r*r' - i*i', r*i' + r'*i) + (r, i) * (r', i') = (r*r' - i*i', r*i' + r'*i) instance Fractional Complex where {-# INLINE (/) #-} - (a, b) / (c, d) - = let den = c^(2 :: Int) + d^(2 :: Int) - re = (a * c + b * d) / den - im = (b * c - a * d) / den - in (re, im) - - fromRational x = (fromRational x, 0) - + (a, b) / (c, d) + = let den = c^(2 :: Int) + d^(2 :: Int) + re = (a * c + b * d) / den + im = (b * c - a * d) / den + in (re, im) + + fromRational x = (fromRational x, 0) + -- | Take the magnitude of a complex number. mag :: Complex -> Double {-# INLINE mag #-} -mag (r, i) = sqrt (r * r + i * i) +mag (r, i) = sqrt (r * r + i * i) -- | Take the argument (phase) of a complex number, in the range [-pi .. pi]. @@ -55,13 +55,13 @@ arg (re, im) = normaliseAngle $ atan2 im re - where normaliseAngle :: Double -> Double - normaliseAngle f - | f < - pi - = normaliseAngle (f + 2 * pi) - - | f > pi - = normaliseAngle (f - 2 * pi) + where normaliseAngle :: Double -> Double + normaliseAngle f + | f < - pi + = normaliseAngle (f + 2 * pi) + + | f > pi + = normaliseAngle (f - 2 * pi) - | otherwise - = f + | otherwise + = f diff --git a/Data/Array/Repa/Algorithms/Convolve.hs b/Data/Array/Repa/Algorithms/Convolve.hs --- a/Data/Array/Repa/Algorithms/Convolve.hs +++ b/Data/Array/Repa/Algorithms/Convolve.hs @@ -12,21 +12,21 @@ -- then use this version instead. -- module Data.Array.Repa.Algorithms.Convolve - ( -- * Arbitrary boundary handling + ( -- * Arbitrary boundary handling convolveP -- * Specialised boundary handling - , GetOut - , outAs - , outClamp - , convolveOutP ) + , GetOut + , outAs + , outClamp + , convolveOutP ) where -import Data.Array.Repa as R +import Data.Array.Repa as R import Data.Array.Repa.Unsafe as R import Data.Array.Repa.Repr.Unboxed as R -import qualified Data.Vector.Unboxed as V -import qualified Data.Array.Repa.Shape as S -import Prelude as P +import qualified Data.Vector.Unboxed as V +import qualified Data.Array.Repa.Shape as S +import Prelude as P -- Plain Convolve ------------------------------------------------------------- @@ -34,70 +34,70 @@ -- which takes a function specifying what value to return when the -- kernel doesn't apply. convolveP - :: (Num a, Unbox a, Monad m) - => (DIM2 -> a) -- ^ Function to get border elements when + :: (Num a, Unbox a, Monad m) + => (DIM2 -> a) -- ^ Function to get border elements when -- the stencil does not apply. - -> Array U DIM2 a -- ^ Stencil to use in the convolution. - -> Array U DIM2 a -- ^ Input image. - -> m (Array U DIM2 a) + -> Array U DIM2 a -- ^ Stencil to use in the convolution. + -> Array U DIM2 a -- ^ Input image. + -> m (Array U DIM2 a) convolveP makeOut kernel image = kernel `deepSeqArray` image `deepSeqArray` computeP $ unsafeTraverse image id update - where + where (Z :. krnHeight :. krnWidth) = extent kernel krnVec = toUnboxed kernel imgSh@(Z :. imgHeight :. imgWidth) = extent image imgVec = toUnboxed image - !krnHeight2 = krnHeight `div` 2 - !krnWidth2 = krnWidth `div` 2 + !krnHeight2 = krnHeight `div` 2 + !krnWidth2 = krnWidth `div` 2 - -- If we're too close to the edge of the input image then - -- we can't apply the stencil because we don't have enough data. - !borderLeft = krnWidth2 - !borderRight = imgWidth - krnWidth2 - 1 - !borderUp = krnHeight2 - !borderDown = imgHeight - krnHeight2 - 1 + -- If we're too close to the edge of the input image then + -- we can't apply the stencil because we don't have enough data. + !borderLeft = krnWidth2 + !borderRight = imgWidth - krnWidth2 - 1 + !borderUp = krnHeight2 + !borderDown = imgHeight - krnHeight2 - 1 - {-# INLINE update #-} - update _ ix@(_ :. j :. i) - | i < borderLeft = makeOut ix - | i > borderRight = makeOut ix - | j < borderUp = makeOut ix - | j > borderDown = makeOut ix - | otherwise = stencil j i + {-# INLINE update #-} + update _ ix@(_ :. j :. i) + | i < borderLeft = makeOut ix + | i > borderRight = makeOut ix + | j < borderUp = makeOut ix + | j > borderDown = makeOut ix + | otherwise = stencil j i - -- The actual stencil function. - {-# INLINE stencil #-} - stencil j i - = let imgStart = S.toIndex imgSh (Z :. j - krnHeight2 :. i - krnWidth2) - in integrate 0 0 0 imgStart 0 + -- The actual stencil function. + {-# INLINE stencil #-} + stencil j i + = let imgStart = S.toIndex imgSh (Z :. j - krnHeight2 :. i - krnWidth2) + in integrate 0 0 0 imgStart 0 - {-# INLINE integrate #-} - integrate !acc !x !y !imgCur !krnCur - | y >= krnHeight - = acc + {-# INLINE integrate #-} + integrate !acc !x !y !imgCur !krnCur + | y >= krnHeight + = acc - | x >= krnWidth - = integrate acc 0 (y + 1) (imgCur + imgWidth - krnWidth) krnCur - - | otherwise - = let imgZ = imgVec `V.unsafeIndex` imgCur - krnZ = krnVec `V.unsafeIndex` krnCur - here = imgZ * krnZ - in integrate (acc + here) (x + 1) y (imgCur + 1) (krnCur + 1) + | x >= krnWidth + = integrate acc 0 (y + 1) (imgCur + imgWidth - krnWidth) krnCur + + | otherwise + = let imgZ = imgVec `V.unsafeIndex` imgCur + krnZ = krnVec `V.unsafeIndex` krnCur + here = imgZ * krnZ + in integrate (acc + here) (x + 1) y (imgCur + 1) (krnCur + 1) {-# INLINE convolveP #-} -- Convolve Out ----------------------------------------------------------------------------------- -- | A function that gets out of range elements from an image. type GetOut a - = (DIM2 -> a) -- ^ The original get function. - -> DIM2 -- ^ The shape of the image. - -> DIM2 -- ^ Index of element we were trying to get. - -> a + = (DIM2 -> a) -- ^ The original get function. + -> DIM2 -- ^ The shape of the image. + -> DIM2 -- ^ Index of element we were trying to get. + -> a -- | Use the provided value for every out-of-range element. @@ -112,70 +112,70 @@ {-# INLINE outClamp #-} outClamp get (_ :. yLen :. xLen) (sh :. j :. i) = clampX j i - where {-# INLINE clampX #-} - clampX !y !x - | x < 0 = clampY y 0 - | x >= xLen = clampY y (xLen - 1) - | otherwise = clampY y x - - {-# INLINE clampY #-} - clampY !y !x - | y < 0 = get (sh :. 0 :. x) - | y >= yLen = get (sh :. (yLen - 1) :. x) - | otherwise = get (sh :. y :. x) + where {-# INLINE clampX #-} + clampX !y !x + | x < 0 = clampY y 0 + | x >= xLen = clampY y (xLen - 1) + | otherwise = clampY y x + + {-# INLINE clampY #-} + clampY !y !x + | y < 0 = get (sh :. 0 :. x) + | y >= yLen = get (sh :. (yLen - 1) :. x) + | otherwise = get (sh :. y :. x) -- | Image-kernel convolution, -- which takes a function specifying what value to use for out-of-range elements. convolveOutP - :: (Num a, Unbox a, Monad m) - => GetOut a -- ^ How to handle out-of-range elements. - -> Array U DIM2 a -- ^ Stencil to use in the convolution. - -> Array U DIM2 a -- ^ Input image. - -> m (Array U DIM2 a) + :: (Num a, Unbox a, Monad m) + => GetOut a -- ^ How to handle out-of-range elements. + -> Array U DIM2 a -- ^ Stencil to use in the convolution. + -> Array U DIM2 a -- ^ Input image. + -> m (Array U DIM2 a) convolveOutP getOut kernel image = kernel `deepSeqArray` image `deepSeqArray` computeP $ unsafeTraverse image id stencil - where + where krnSh@(Z :. krnHeight :. krnWidth) = extent kernel imgSh@(Z :. imgHeight :. imgWidth) = extent image - !krnHeight2 = krnHeight `div` 2 - !krnWidth2 = krnWidth `div` 2 - !krnSize = S.size krnSh + !krnHeight2 = krnHeight `div` 2 + !krnWidth2 = krnWidth `div` 2 + !krnSize = S.size krnSh - -- If we're too close to the edge of the input image then - -- we can't apply the stencil because we don't have enough data. - !borderLeft = krnWidth2 - !borderRight = imgWidth - krnWidth2 - 1 - !borderUp = krnHeight2 - !borderDown = imgHeight - krnHeight2 - 1 + -- If we're too close to the edge of the input image then + -- we can't apply the stencil because we don't have enough data. + !borderLeft = krnWidth2 + !borderRight = imgWidth - krnWidth2 - 1 + !borderUp = krnHeight2 + !borderDown = imgHeight - krnHeight2 - 1 - -- The actual stencil function. - {-# INLINE stencil #-} - stencil get (_ :. j :. i) - = let - {-# INLINE get' #-} - get' ix@(_ :. y :. x) - | x < borderLeft = getOut get imgSh ix - | x > borderRight = getOut get imgSh ix - | y < borderUp = getOut get imgSh ix - | y > borderDown = getOut get imgSh ix - | otherwise = get ix + -- The actual stencil function. + {-# INLINE stencil #-} + stencil get (_ :. j :. i) + = let + {-# INLINE get' #-} + get' ix@(_ :. y :. x) + | x < borderLeft = getOut get imgSh ix + | x > borderRight = getOut get imgSh ix + | y < borderUp = getOut get imgSh ix + | y > borderDown = getOut get imgSh ix + | otherwise = get ix - !ikrnWidth' = i - krnWidth2 - !jkrnHeight' = j - krnHeight2 + !ikrnWidth' = i - krnWidth2 + !jkrnHeight' = j - krnHeight2 - {-# INLINE integrate #-} - integrate !count !acc - | count == krnSize = acc - | otherwise - = let !ix@(sh :. y :. x) = S.fromIndex krnSh count - !ix' = sh :. y + jkrnHeight' :. x + ikrnWidth' - !here = kernel `unsafeIndex` ix * (get' ix') - in integrate (count + 1) (acc + here) + {-# INLINE integrate #-} + integrate !count !acc + | count == krnSize = acc + | otherwise + = let !ix@(sh :. y :. x) = S.fromIndex krnSh count + !ix' = sh :. y + jkrnHeight' :. x + ikrnWidth' + !here = kernel `unsafeIndex` ix * (get' ix') + in integrate (count + 1) (acc + here) - in integrate 0 0 + in integrate 0 0 {-# INLINE convolveOutP #-} diff --git a/Data/Array/Repa/Algorithms/DFT.hs b/Data/Array/Repa/Algorithms/DFT.hs --- a/Data/Array/Repa/Algorithms/DFT.hs +++ b/Data/Array/Repa/Algorithms/DFT.hs @@ -14,89 +14,89 @@ -- -- You can also compute single values of the transform using `dftWithRootsSingle`. module Data.Array.Repa.Algorithms.DFT - ( dftP - , idftP - , dftWithRootsP - , dftWithRootsSingleS) + ( dftP + , idftP + , dftWithRootsP + , dftWithRootsSingleS) where -import Data.Array.Repa.Algorithms.DFT.Roots -import Data.Array.Repa.Algorithms.Complex -import Data.Array.Repa as R -import Prelude as P +import Data.Array.Repa.Algorithms.DFT.Roots as R +import Data.Array.Repa.Algorithms.Complex as R +import Data.Array.Repa as R +import Prelude as P -- | Compute the DFT along the low order dimension of an array. -dftP :: (Shape sh, Monad m) - => Array U (sh :. Int) Complex - -> m (Array U (sh :. Int) Complex) +dftP :: (Shape sh, Monad m) + => Array U (sh :. Int) Complex + -> m (Array U (sh :. Int) Complex) dftP v - = do rofu <- calcRootsOfUnityP (extent v) + = do rofu <- calcRootsOfUnityP (extent v) dftWithRootsP rofu v {-# INLINE dftP #-} -- | Compute the inverse DFT along the low order dimension of an array. -idftP :: (Shape sh, Monad m) - => Array U (sh :. Int) Complex - -> m (Array U (sh :. Int) Complex) +idftP :: (Shape sh, Monad m) + => Array U (sh :. Int) Complex + -> m (Array U (sh :. Int) Complex) idftP v - = do let _ :. len = extent v - let scale = (fromIntegral len, 0) - rofu <- calcInverseRootsOfUnityP (extent v) + = do let _ :. len = extent v + let scale = (fromIntegral len, 0) + rofu <- calcInverseRootsOfUnityP (extent v) roots <- dftWithRootsP rofu v computeP $ R.map (/ scale) roots {-# INLINE idftP #-} -- | Generic function for computation of forward or inverse DFT. --- This function is also useful if you transform many arrays with the same extent, --- and don't want to recompute the roots for each one. --- The extent of the given roots must match that of the input array, else `error`. +-- This function is also useful if you transform many arrays with the same extent, +-- and don't want to recompute the roots for each one. +-- The extent of the given roots must match that of the input array, else `error`. dftWithRootsP - :: (Shape sh, Monad m) - => Array U (sh :. Int) Complex -- ^ Roots of unity. - -> Array U (sh :. Int) Complex -- ^ Input array. - -> m (Array U (sh :. Int) Complex) + :: (Shape sh, Monad m) + => Array U (sh :. Int) Complex -- ^ Roots of unity. + -> Array U (sh :. Int) Complex -- ^ Input array. + -> m (Array U (sh :. Int) Complex) dftWithRootsP rofu arr - | _ :. rLen <- extent rofu - , _ :. vLen <- extent arr - , rLen /= vLen - = error $ "dftWithRoots: length of vector (" P.++ show vLen P.++ ")" - P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")" + | _ :. rLen <- extent rofu + , _ :. vLen <- extent arr + , rLen /= vLen + = error $ "dftWithRoots: length of vector (" P.++ show vLen P.++ ")" + P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")" - | otherwise - = computeP $ traverse arr id (\_ k -> dftWithRootsSingleS rofu arr k) -{-# INLINE dftWithRootsP #-} + | otherwise + = computeP $ R.traverse arr id (\_ k -> dftWithRootsSingleS rofu arr k) +{-# INLINE dftWithRootsP #-} -- | Compute a single value of the DFT. --- The extent of the given roots must match that of the input array, else `error`. +-- The extent of the given roots must match that of the input array, else `error`. dftWithRootsSingleS - :: Shape sh - => Array U (sh :. Int) Complex -- ^ Roots of unity. - -> Array U (sh :. Int) Complex -- ^ Input array. - -> (sh :. Int) -- ^ Index of the value we want. - -> Complex + :: Shape sh + => Array U (sh :. Int) Complex -- ^ Roots of unity. + -> Array U (sh :. Int) Complex -- ^ Input array. + -> (sh :. Int) -- ^ Index of the value we want. + -> Complex dftWithRootsSingleS rofu arrX (_ :. k) - | _ :. rLen <- extent rofu - , _ :. vLen <- extent arrX - , rLen /= vLen - = error $ "dftWithRootsSingle: length of vector (" P.++ show vLen P.++ ")" - P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")" + | _ :. rLen <- extent rofu + , _ :. vLen <- extent arrX + , rLen /= vLen + = error $ "dftWithRootsSingle: length of vector (" P.++ show vLen P.++ ")" + P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")" - | otherwise - = let sh@(_ :. len) = extent arrX + | otherwise + = let sh@(_ :. len) = extent arrX - -- All the roots we need to multiply with. - wroots = fromFunction sh elemFn - elemFn (sh' :. n) - = rofu ! (sh' :. (k * n) `mod` len) + -- All the roots we need to multiply with. + wroots = fromFunction sh elemFn + elemFn (sh' :. n) + = rofu ! (sh' :. (k * n) `mod` len) - in R.sumAllS $ R.zipWith (*) arrX wroots + in R.sumAllS $ R.zipWith (*) arrX wroots {-# INLINE dftWithRootsSingleS #-} diff --git a/Data/Array/Repa/Algorithms/DFT/Center.hs b/Data/Array/Repa/Algorithms/DFT/Center.hs --- a/Data/Array/Repa/Algorithms/DFT/Center.hs +++ b/Data/Array/Repa/Algorithms/DFT/Center.hs @@ -2,12 +2,12 @@ -- | Applying these transforms to the input of a DFT causes the output -- to be centered so that the zero frequency is in the middle. module Data.Array.Repa.Algorithms.DFT.Center - ( center1d - , center2d - , center3d) + ( center1d + , center2d + , center3d) where -import Data.Array.Repa -import Data.Array.Repa.Algorithms.Complex +import Data.Array.Repa as R +import Data.Array.Repa.Algorithms.Complex as R -- | Apply the centering transform to a vector. center1d @@ -15,8 +15,8 @@ => Array r DIM1 Complex -> Array D DIM1 Complex {-# INLINE center1d #-} center1d arr - = traverse arr id - (\get ix@(_ :. x) -> ((-1) ^ x) * get ix) + = R.traverse arr id + (\get ix@(_ :. x) -> ((-1) ^ x) * get ix) -- | Apply the centering transform to a matrix. @@ -25,8 +25,8 @@ => Array r DIM2 Complex -> Array D DIM2 Complex {-# INLINE center2d #-} center2d arr - = traverse arr id - (\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix) + = R.traverse arr id + (\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix) -- | Apply the centering transform to a 3d array. @@ -35,5 +35,5 @@ => Array r DIM3 Complex -> Array D DIM3 Complex {-# INLINE center3d #-} center3d arr - = traverse arr id - (\get ix@(_ :. z :. y :. x) -> ((-1) ^ (z + y + x)) * get ix) + = R.traverse arr id + (\get ix@(_ :. z :. y :. x) -> ((-1) ^ (z + y + x)) * get ix) diff --git a/Data/Array/Repa/Algorithms/DFT/Roots.hs b/Data/Array/Repa/Algorithms/DFT/Roots.hs --- a/Data/Array/Repa/Algorithms/DFT/Roots.hs +++ b/Data/Array/Repa/Algorithms/DFT/Roots.hs @@ -2,8 +2,8 @@ -- | Calculation of roots of unity for the forward and inverse DFT\/FFT. module Data.Array.Repa.Algorithms.DFT.Roots - ( calcRootsOfUnityP - , calcInverseRootsOfUnityP) + ( calcRootsOfUnityP + , calcInverseRootsOfUnityP) where import Data.Array.Repa import Data.Array.Repa.Algorithms.Complex @@ -11,33 +11,33 @@ -- | Calculate roots of unity for the forward transform. calcRootsOfUnityP - :: (Shape sh, Monad m) - => (sh :. Int) -- ^ Length of lowest dimension of result. - -> m (Array U (sh :. Int) Complex) + :: (Shape sh, Monad m) + => (sh :. Int) -- ^ Length of lowest dimension of result. + -> m (Array U (sh :. Int) Complex) calcRootsOfUnityP sh@(_ :. n) = computeP $ fromFunction sh f where f :: Shape sh => (sh :. Int) -> Complex f (_ :. i) - = ( cos (2 * pi * (fromIntegral i) / len) - , - sin (2 * pi * (fromIntegral i) / len)) + = ( cos (2 * pi * (fromIntegral i) / len) + , - sin (2 * pi * (fromIntegral i) / len)) - len = fromIntegral n + len = fromIntegral n -- | Calculate roots of unity for the inverse transform. calcInverseRootsOfUnityP - :: (Shape sh, Monad m) - => (sh :. Int) -- ^ Length of lowest dimension of result. - -> m (Array U (sh :. Int) Complex) + :: (Shape sh, Monad m) + => (sh :. Int) -- ^ Length of lowest dimension of result. + -> m (Array U (sh :. Int) Complex) calcInverseRootsOfUnityP sh@(_ :. n) = computeP $ fromFunction sh f where f :: Shape sh => (sh :. Int) -> Complex f (_ :. i) - = ( cos (2 * pi * (fromIntegral i) / len) - , sin (2 * pi * (fromIntegral i) / len)) + = ( cos (2 * pi * (fromIntegral i) / len) + , sin (2 * pi * (fromIntegral i) / len)) - len = fromIntegral n + len = fromIntegral n diff --git a/Data/Array/Repa/Algorithms/FFT.hs b/Data/Array/Repa/Algorithms/FFT.hs --- a/Data/Array/Repa/Algorithms/FFT.hs +++ b/Data/Array/Repa/Algorithms/FFT.hs @@ -8,82 +8,82 @@ -- 50x slower than FFTW in estimate mode. -- module Data.Array.Repa.Algorithms.FFT - ( Mode(..) - , isPowerOfTwo - , fft3dP - , fft2dP - , fft1dP) + ( Mode(..) + , isPowerOfTwo + , fft3dP + , fft2dP + , fft1dP) where import Data.Array.Repa.Algorithms.Complex -import Data.Array.Repa as R +import Data.Array.Repa as R import Data.Array.Repa.Eval as R import Data.Array.Repa.Unsafe as R import Prelude as P data Mode - = Forward - | Reverse - | Inverse - deriving (Show, Eq) + = Forward + | Reverse + | Inverse + deriving (Show, Eq) signOfMode :: Mode -> Double signOfMode mode = case mode of - Forward -> (-1) - Reverse -> 1 - Inverse -> 1 + Forward -> (-1) + Reverse -> 1 + Inverse -> 1 {-# INLINE signOfMode #-} -- | Check if an `Int` is a power of two. isPowerOfTwo :: Int -> Bool isPowerOfTwo n - | 0 <- n = True - | 2 <- n = True - | n `mod` 2 == 0 = isPowerOfTwo (n `div` 2) - | otherwise = False + | 0 <- n = True + | 2 <- n = True + | n `mod` 2 == 0 = isPowerOfTwo (n `div` 2) + | otherwise = False {-# INLINE isPowerOfTwo #-} -- 3D Transform ----------------------------------------------------------------------------------- -- | Compute the DFT of a 3d array. Array dimensions must be powers of two else `error`. -fft3dP :: (Source r Complex, Monad m) +fft3dP :: (Source r Complex, Monad m) => Mode - -> Array r DIM3 Complex - -> m (Array U DIM3 Complex) + -> Array r DIM3 Complex + -> m (Array U DIM3 Complex) fft3dP mode arr - = let _ :. depth :. height :. width = extent arr - !sign = signOfMode mode - !scale = fromIntegral (depth * width * height) - - in if not (isPowerOfTwo depth && isPowerOfTwo height && isPowerOfTwo width) - then error $ unlines - [ "Data.Array.Repa.Algorithms.FFT: fft3d" - , " Array dimensions must be powers of two," - , " but the provided array is " - P.++ show height P.++ "x" P.++ show width P.++ "x" P.++ show depth ] - - else arr `deepSeqArray` - case mode of - Forward -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr - Reverse -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr - Inverse -> computeP - $ R.map (/ scale) - $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr + = let _ :. depth :. height :. width = extent arr + !sign = signOfMode mode + !scale = fromIntegral (depth * width * height) + + in if not (isPowerOfTwo depth && isPowerOfTwo height && isPowerOfTwo width) + then error $ unlines + [ "Data.Array.Repa.Algorithms.FFT: fft3d" + , " Array dimensions must be powers of two," + , " but the provided array is " + P.++ show height P.++ "x" P.++ show width P.++ "x" P.++ show depth ] + + else arr `deepSeqArray` + case mode of + Forward -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr + Reverse -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr + Inverse -> computeP + $ R.map (/ scale) + $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr {-# INLINE fft3dP #-} fftTrans3d - :: Source r Complex - => Double - -> Array r DIM3 Complex - -> Array U DIM3 Complex + :: Source r Complex + => Double + -> Array r DIM3 Complex + -> Array U DIM3 Complex fftTrans3d sign arr - = let (sh :. len) = extent arr - in suspendedComputeP $ rotate3d $ fft sign sh len arr + = let (sh :. len) = extent arr + in suspendedComputeP $ rotate3d $ fft sign sh len arr {-# INLINE fftTrans3d #-} @@ -92,83 +92,83 @@ => Array r DIM3 Complex -> Array D DIM3 Complex rotate3d arr = backpermute (sh :. m :. k :. l) f arr - where (sh :. k :. l :. m) = extent arr - f (sh' :. m' :. k' :. l') = sh' :. k' :. l' :. m' + where (sh :. k :. l :. m) = extent arr + f (sh' :. m' :. k' :. l') = sh' :. k' :. l' :. m' {-# INLINE rotate3d #-} -- Matrix Transform ------------------------------------------------------------------------------- -- | Compute the DFT of a matrix. Array dimensions must be powers of two else `error`. -fft2dP :: (Source r Complex, Monad m) +fft2dP :: (Source r Complex, Monad m) => Mode - -> Array r DIM2 Complex - -> m (Array U DIM2 Complex) + -> Array r DIM2 Complex + -> m (Array U DIM2 Complex) fft2dP mode arr - = let _ :. height :. width = extent arr - sign = signOfMode mode - scale = fromIntegral (width * height) - - in if not (isPowerOfTwo height && isPowerOfTwo width) - then error $ unlines - [ "Data.Array.Repa.Algorithms.FFT: fft2d" - , " Array dimensions must be powers of two," - , " but the provided array is " P.++ show height P.++ "x" P.++ show width ] - - else arr `deepSeqArray` - case mode of - Forward -> now $ fftTrans2d sign $ fftTrans2d sign arr - Reverse -> now $ fftTrans2d sign $ fftTrans2d sign arr - Inverse -> computeP $ R.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr + = let _ :. height :. width = extent arr + sign = signOfMode mode + scale = fromIntegral (width * height) + + in if not (isPowerOfTwo height && isPowerOfTwo width) + then error $ unlines + [ "Data.Array.Repa.Algorithms.FFT: fft2d" + , " Array dimensions must be powers of two," + , " but the provided array is " P.++ show height P.++ "x" P.++ show width ] + + else arr `deepSeqArray` + case mode of + Forward -> now $ fftTrans2d sign $ fftTrans2d sign arr + Reverse -> now $ fftTrans2d sign $ fftTrans2d sign arr + Inverse -> computeP $ R.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr {-# INLINE fft2dP #-} fftTrans2d - :: Source r Complex - => Double - -> Array r DIM2 Complex - -> Array U DIM2 Complex + :: Source r Complex + => Double + -> Array r DIM2 Complex + -> Array U DIM2 Complex fftTrans2d sign arr - = let (sh :. len) = extent arr - in suspendedComputeP $ transpose $ fft sign sh len arr + = let (sh :. len) = extent arr + in suspendedComputeP $ transpose $ fft sign sh len arr {-# INLINE fftTrans2d #-} -- Vector Transform ------------------------------------------------------------------------------- -- | Compute the DFT of a vector. Array dimensions must be powers of two else `error`. -fft1dP :: (Source r Complex, Monad m) +fft1dP :: (Source r Complex, Monad m) => Mode - -> Array r DIM1 Complex - -> m (Array U DIM1 Complex) + -> Array r DIM1 Complex + -> m (Array U DIM1 Complex) fft1dP mode arr - = let _ :. len = extent arr - sign = signOfMode mode - scale = fromIntegral len - - in if not $ isPowerOfTwo len - then error $ unlines + = let _ :. len = extent arr + sign = signOfMode mode + scale = fromIntegral len + + in if not $ isPowerOfTwo len + then error $ unlines [ "Data.Array.Repa.Algorithms.FFT: fft1d" - , " Array dimensions must be powers of two, " + , " Array dimensions must be powers of two, " , " but the provided array is " P.++ show len ] - - else arr `deepSeqArray` - case mode of - Forward -> now $ fftTrans1d sign arr - Reverse -> now $ fftTrans1d sign arr - Inverse -> computeP $ R.map (/ scale) $ fftTrans1d sign arr + + else arr `deepSeqArray` + case mode of + Forward -> now $ fftTrans1d sign arr + Reverse -> now $ fftTrans1d sign arr + Inverse -> computeP $ R.map (/ scale) $ fftTrans1d sign arr {-# INLINE fft1dP #-} fftTrans1d - :: Source r Complex - => Double - -> Array r DIM1 Complex - -> Array U DIM1 Complex + :: Source r Complex + => Double + -> Array r DIM1 Complex + -> Array U DIM1 Complex fftTrans1d sign arr - = let (sh :. len) = extent arr - in fft sign sh len arr + = let (sh :. len) = extent arr + in fft sign sh len arr {-# INLINE fftTrans1d #-} @@ -180,37 +180,37 @@ fft !sign !sh !lenVec !vec = go lenVec 0 1 - where go !len !offset !stride - | len == 2 - = suspendedComputeP $ fromFunction (sh :. 2) swivel - - | otherwise - = combine len - (go (len `div` 2) offset (stride * 2)) - (go (len `div` 2) (offset + stride) (stride * 2)) + where go !len !offset !stride + | len == 2 + = suspendedComputeP $ fromFunction (sh :. 2) swivel + + | otherwise + = combine len + (go (len `div` 2) offset (stride * 2)) + (go (len `div` 2) (offset + stride) (stride * 2)) - where swivel (sh' :. ix) - = case ix of - 0 -> (vec `unsafeIndex` (sh' :. offset)) + (vec `unsafeIndex` (sh' :. (offset + stride))) - 1 -> (vec `unsafeIndex` (sh' :. offset)) - (vec `unsafeIndex` (sh' :. (offset + stride))) + where swivel (sh' :. ix) + = case ix of + 0 -> (vec `unsafeIndex` (sh' :. offset)) + (vec `unsafeIndex` (sh' :. (offset + stride))) + 1 -> (vec `unsafeIndex` (sh' :. offset)) - (vec `unsafeIndex` (sh' :. (offset + stride))) - {-# INLINE combine #-} - combine !len' evens odds - = evens `deepSeqArray` odds `deepSeqArray` - let odds' = unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) - in suspendedComputeP $ (evens +^ odds') R.++ (evens -^ odds') + {-# INLINE combine #-} + combine !len' evens odds + = evens `deepSeqArray` odds `deepSeqArray` + let odds' = unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) + in suspendedComputeP $ (evens +^ odds') R.++ (evens -^ odds') {-# INLINE fft #-} -- Compute a twiddle factor. twiddle :: Double - -> Int -- index - -> Int -- length - -> Complex + -> Int -- index + -> Int -- length + -> Complex twiddle sign k' n' - = (cos (2 * pi * k / n), sign * sin (2 * pi * k / n)) - where k = fromIntegral k' - n = fromIntegral n' + = (cos (2 * pi * k / n), sign * sin (2 * pi * k / n)) + where k = fromIntegral k' + n = fromIntegral n' {-# INLINE twiddle #-} diff --git a/Data/Array/Repa/Algorithms/Matrix.hs b/Data/Array/Repa/Algorithms/Matrix.hs --- a/Data/Array/Repa/Algorithms/Matrix.hs +++ b/Data/Array/Repa/Algorithms/Matrix.hs @@ -15,7 +15,7 @@ , col -- * Matrix Multiplication. - , mmultP, mmultS + , mmultP, mmultS -- * Transposition. , transpose2P, transpose2S diff --git a/Data/Array/Repa/Algorithms/Randomish.hs b/Data/Array/Repa/Algorithms/Randomish.hs --- a/Data/Array/Repa/Algorithms/Randomish.hs +++ b/Data/Array/Repa/Algorithms/Randomish.hs @@ -1,17 +1,17 @@ {-# LANGUAGE BangPatterns #-} module Data.Array.Repa.Algorithms.Randomish - ( randomishIntArray - , randomishIntVector - , randomishDoubleArray - , randomishDoubleVector) + ( randomishIntArray + , randomishIntVector + , randomishDoubleArray + , randomishDoubleVector) where import Data.Word -import Data.Vector.Unboxed (Vector) -import Data.Array.Repa as R -import qualified Data.Vector.Unboxed.Mutable as MV -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Generic as G +import Data.Vector.Unboxed (Vector) +import Data.Array.Repa as R +import qualified Data.Vector.Unboxed.Mutable as MV +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Generic as G -- | Use the ''minimal standard'' Lehmer generator to quickly generate some random @@ -27,89 +27,89 @@ -- Communications of the ACM, Oct 1988, Volume 31, Number 10. -- randomishIntArray - :: Shape sh - => sh -- ^ Shape of array - -> Int -- ^ Minumum value in output. - -> Int -- ^ Maximum value in output. - -> Int -- ^ Random seed. - -> Array U sh Int -- ^ Array of randomish numbers. + :: Shape sh + => sh -- ^ Shape of array + -> Int -- ^ Minumum value in output. + -> Int -- ^ Maximum value in output. + -> Int -- ^ Random seed. + -> Array U sh Int -- ^ Array of randomish numbers. randomishIntArray !sh !valMin !valMax !seed - = fromUnboxed sh $ randomishIntVector (R.size sh) valMin valMax seed + = fromUnboxed sh $ randomishIntVector (R.size sh) valMin valMax seed randomishIntVector - :: Int -- ^ Length of vector. - -> Int -- ^ Minumum value in output. - -> Int -- ^ Maximum value in output. - -> Int -- ^ Random seed. - -> Vector Int -- ^ Vector of randomish numbers. + :: Int -- ^ Length of vector. + -> Int -- ^ Minumum value in output. + -> Int -- ^ Maximum value in output. + -> Int -- ^ Random seed. + -> Vector Int -- ^ Vector of randomish numbers. randomishIntVector !len !valMin' !valMax' !seed' - = let -- a magic number - -- (don't change it, the randomness depends on this specific number). - multiplier :: Word64 - multiplier = 16807 + = let -- a magic number + -- (don't change it, the randomness depends on this specific number). + multiplier :: Word64 + multiplier = 16807 - -- a merzenne prime - -- (don't change it, the randomness depends on this specific number). - modulus :: Word64 - modulus = 2^(31 :: Integer) - 1 + -- a merzenne prime + -- (don't change it, the randomness depends on this specific number). + modulus :: Word64 + modulus = 2^(31 :: Integer) - 1 - -- if the seed is 0 all the numbers in the sequence are also 0. - seed - | seed' == 0 = 1 - | otherwise = seed' + -- if the seed is 0 all the numbers in the sequence are also 0. + seed + | seed' == 0 = 1 + | otherwise = seed' - !valMin = fromIntegral valMin' - !valMax = fromIntegral valMax' + 1 - !range = valMax - valMin + !valMin = fromIntegral valMin' + !valMax = fromIntegral valMax' + 1 + !range = valMax - valMin - {-# INLINE f #-} - f x = multiplier * x `mod` modulus + {-# INLINE f #-} + f x = multiplier * x `mod` modulus in G.create - $ do - vec <- MV.new len + $ do + vec <- MV.new len - let go !ix !x - | ix == len = return () - | otherwise - = do let x' = f x - MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin - go (ix + 1) x' + let go !ix !x + | ix == len = return () + | otherwise + = do let x' = f x + MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin + go (ix + 1) x' - go 0 (f $ f $ f $ fromIntegral seed) - return vec + go 0 (f $ f $ f $ fromIntegral seed) + return vec -- | Generate some randomish doubles with terrible statistical properties. -- This just takes randomish ints then scales them, so there's not much randomness in low-order bits. randomishDoubleArray - :: Shape sh - => sh -- ^ Shape of array - -> Double -- ^ Minumum value in output. - -> Double -- ^ Maximum value in output. - -> Int -- ^ Random seed. - -> Array U sh Double -- ^ Array of randomish numbers. + :: Shape sh + => sh -- ^ Shape of array + -> Double -- ^ Minumum value in output. + -> Double -- ^ Maximum value in output. + -> Int -- ^ Random seed. + -> Array U sh Double -- ^ Array of randomish numbers. randomishDoubleArray !sh !valMin !valMax !seed - = fromUnboxed sh $ randomishDoubleVector (R.size sh) valMin valMax seed + = fromUnboxed sh $ randomishDoubleVector (R.size sh) valMin valMax seed -- | Generate some randomish doubles with terrible statistical properties. -- This just takes randmish ints then scales them, so there's not much randomness in low-order bits. randomishDoubleVector - :: Int -- ^ Length of vector - -> Double -- ^ Minimum value in output - -> Double -- ^ Maximum value in output - -> Int -- ^ Random seed. - -> Vector Double -- ^ Vector of randomish doubles. + :: Int -- ^ Length of vector + -> Double -- ^ Minimum value in output + -> Double -- ^ Maximum value in output + -> Int -- ^ Random seed. + -> Vector Double -- ^ Vector of randomish doubles. randomishDoubleVector !len !valMin !valMax !seed - = let range = valMax - valMin + = let range = valMax - valMin - mx = 2^(30 :: Integer) - 1 - mxf = fromIntegral (mx :: Integer) - ints = randomishIntVector len 0 mx seed - - in V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints + mx = 2^(30 :: Integer) - 1 + mxf = fromIntegral (mx :: Integer) + ints = randomishIntVector len 0 mx seed + + in V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints diff --git a/repa-algorithms.cabal b/repa-algorithms.cabal --- a/repa-algorithms.cabal +++ b/repa-algorithms.cabal @@ -1,5 +1,5 @@ Name: repa-algorithms -Version: 3.3.1.2 +Version: 3.4.0.1 License: BSD3 License-file: LICENSE Author: The DPH Team @@ -18,9 +18,9 @@ Library Build-Depends: - base == 4.7.*, + base == 4.8.*, vector == 0.10.*, - repa == 3.3.1.* + repa == 3.4.0.* ghc-options: -Wall -fno-warn-missing-signatures @@ -44,7 +44,6 @@ StandaloneDeriving ScopedTypeVariables PatternGuards - OverlappingInstances Exposed-modules: Data.Array.Repa.Algorithms.DFT.Center