repa-algorithms 3.3.1.2 → 3.4.0.1

raw patch · 10 files changed

+425/−426 lines, 10 filesdep ~basedep ~repaPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base, repa

API changes (from Hackage documentation)

Files

• Data/Array/Repa/Algorithms/ColorRamp.hs +37/−37
• Data/Array/Repa/Algorithms/Complex.hs +28/−28
• Data/Array/Repa/Algorithms/Convolve.hs +99/−99
• Data/Array/Repa/Algorithms/DFT.hs +51/−51
• Data/Array/Repa/Algorithms/DFT/Center.hs +11/−11
• Data/Array/Repa/Algorithms/DFT/Roots.hs +14/−14
• Data/Array/Repa/Algorithms/FFT.hs +115/−115
• Data/Array/Repa/Algorithms/Matrix.hs +1/−1
• Data/Array/Repa/Algorithms/Randomish.hs +66/−66
• repa-algorithms.cabal +3/−4

@@ -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 

@@ -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

@@ -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 #-} 

@@ -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 #-}  

@@ -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)

@@ -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

@@ -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 #-} 

@@ -15,7 +15,7 @@         , col            -- * Matrix Multiplication.-	, mmultP,      mmultS+        , mmultP,      mmultS            -- * Transposition.         , transpose2P, transpose2S

@@ -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

@@ -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