repa-algorithms 1.1.0.0 → 2.0.0.1

raw patch · 10 files changed

+567/−231 lines, 10 filesdep +vectordep −dph-basedep ~repaPVP ok

version bump matches the API change (PVP)

Dependencies added: vector

Dependencies removed: dph-base

Dependency ranges changed: repa

API changes (from Hackage documentation)

Files

• Data/Array/Repa/Algorithms/Complex.hs +31/−17
• Data/Array/Repa/Algorithms/Convolve.hs +167/−0
• Data/Array/Repa/Algorithms/DFT.hs +7/−6
• Data/Array/Repa/Algorithms/DFT/Center.hs +17/−15
• Data/Array/Repa/Algorithms/DFT/Roots.hs +6/−4
• Data/Array/Repa/Algorithms/FFT.hs +144/−171
• Data/Array/Repa/Algorithms/Iterate.hs +45/−0
• Data/Array/Repa/Algorithms/Matrix.hs +22/−11
• Data/Array/Repa/Algorithms/Randomish.hs +115/−0
• repa-algorithms.cabal +13/−7

@@ -1,44 +1,58 @@-{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}+{-# LANGUAGE TypeOperators, TypeSynonymInstances, FlexibleInstances #-}  -- | Strict complex doubles. module Data.Array.Repa.Algorithms.Complex 	( Complex 	, mag-	, arg-	, (:*:)(..))+	, arg) where-import 	Data.Array.Parallel.Base ((:*:)(..)) --- | Strict complex doubles.++-- | Complex doubles. type Complex -	= Double :*: Double+	= (Double, Double)  instance Num Complex where-  abs x				= (mag x) :*: 0-  signum (re :*: _)		= signum re :*: 0-  fromInteger n			= fromInteger n :*: 0.0-  (r :*: i) + (r' :*: i')	= r+r' :*: i+i'-  (r :*: i) - (r' :*: i')	= r-r' :*: i-i'-  (r :*: i) * (r' :*: i')	= r*r' - i*i' :*: r*i' + r'*i +  {-# INLINE abs #-}+  abs x			= (mag x, 0) +  {-# INLINE signum #-}+  signum (re, _)	= (signum re, 0)++  {-# INLINE fromInteger #-}+  fromInteger n		= (fromInteger n, 0.0)++  {-# INLINE (+) #-}+  (r, i) + (r', i')	= (r+r', i+i')++  {-# INLINE (-) #-}+  (r, i) - (r', i')	= (r-r', i-i')++  {-# INLINE (*) #-}+  (r, i) * (r', i')	= (r*r' - i*i', r*i' + r'*i)++ instance Fractional Complex where-  (a :*: b) / (c :*: d)		+  {-# 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+	  in	(re, im) 	-  fromRational x	= fromRational x :*: 0+  fromRational x	= (fromRational x, 0) 	 -- | Take the magnitude of a complex number. mag :: Complex -> Double-mag (r :*: i)	= sqrt (r * r + i * i)+{-# INLINE mag #-}+mag (r, i)	= sqrt (r * r + i * i)   -- | Take the argument (phase) of a complex number, in the range [-pi .. pi]. arg :: Complex -> Double-arg (re :*: im)+{-# INLINE arg #-}+arg (re, im)  = normaliseAngle $ atan2 im re   where 	normaliseAngle :: Double -> Double

@@ -0,0 +1,167 @@+{-# LANGUAGE BangPatterns, PackageImports #-}+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}++-- | Old support for stencil based convolutions. +--+--   NOTE: This is slated to be merged with the new Stencil support in the next version+--         of Repa. We'll still expose the `convolve` function though.+--+module Data.Array.Repa.Algorithms.Convolve+	( convolve++	, GetOut+	, outAs+	, outClamp+	, convolveOut )+where+import Data.Array.Repa 					as A+import qualified Data.Vector.Unboxed			as V+import qualified Data.Array.Repa.Shape			as S+import Prelude						as P+++-- Plain Convolve ---------------------------------------------------------------------------------+-- | Image-kernel convolution,+--   which takes a function specifying what value to return when the kernel doesn't apply.+convolve+	:: (Elt a, Num a)+	=> (DIM2 -> a) 		-- ^ Use this function to get border elements where the kernel apply.+	-> Array DIM2 a		-- ^ Kernel to use in the convolution.+	-> Array DIM2 a		-- ^ Input image.+	-> Array DIM2 a++{-# INLINE convolve #-}+convolve makeOut+ 	kernel@(Array       (_ :. krnHeight :. krnWidth) [Region RangeAll (GenManifest krnVec)])+  	 image@(Array imgSh@(_ :. imgHeight :. imgWidth) [Region RangeAll (GenManifest imgVec)])++ = kernel `deepSeqArray` image `deepSeqArray` +   force $ unsafeTraverse image id update+ where	+	!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++	{-# 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++	{-# 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)+++-- 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+++-- | Use the provided value for every out-of-range element.+outAs :: a -> GetOut a+{-# INLINE outAs #-}+outAs x _ _ _ = x+++-- | If the requested element is out of range use+--   the closest one from the real image.+outClamp :: GetOut a+{-# 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)+++-- | Image-kernel convolution, +--   which takes a function specifying what value to use for out-of-range elements.+convolveOut+	:: (Elt a, Num a)+	=> GetOut a		-- ^ Use this fn to get out of range elements.+	-> Array DIM2 a		-- ^ Kernel+	-> Array DIM2 a		-- ^ Image+	-> Array DIM2 a++{-# INLINE convolveOut #-}+convolveOut getOut+ 	kernel@(Array krnSh@(_ :. krnHeight :. krnWidth) _)+  	 image@(Array imgSh@(_ :. imgHeight :. imgWidth) _)++ = kernel `deepSeqArray` image `deepSeqArray` +   force $ unsafeTraverse image id stencil+ where	+	!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++	-- 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++		{-# 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+

@@ -22,6 +22,7 @@ import Data.Array.Repa.Algorithms.DFT.Roots import Data.Array.Repa.Algorithms.Complex import Data.Array.Repa				as A+import Prelude					as P  -- | Compute the DFT along the low order dimension of an array. dft 	:: forall sh@@ -42,7 +43,7 @@  idft v  = let	_ :. len	= extent v-	scale		= fromIntegral len :*: 0+	scale		= (fromIntegral len, 0) 	rofu		= calcInverseRootsOfUnity (extent v)    in	force $ A.map (/ scale) $ dftWithRoots rofu v @@ -62,8 +63,8 @@ 	| _ :. rLen 	<- extent rofu 	, _ :. vLen 	<- extent arr 	, rLen /= vLen-	= error $  "dftWithRoots: length of vector (" ++ show vLen ++ ")"-		++ " does not match the length of the roots (" ++ show rLen ++ ")"+	= error $    "dftWithRoots: length of vector (" P.++ show vLen P.++ ")"+		P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")"  	| otherwise 	= traverse arr id (\_ k -> dftWithRootsSingle rofu arr k)@@ -84,8 +85,8 @@ 	| _ :. rLen 	<- extent rofu 	, _ :. vLen 	<- extent arrX 	, rLen /= vLen-	= error $  "dftWithRootsSingle: length of vector (" ++ show vLen ++ ")"-		++ " does not match the length of the roots (" ++ show rLen ++ ")"+	= 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@@ -93,7 +94,7 @@ 		-- All the roots we need to multiply with. 		wroots		= fromFunction sh elemFn 		elemFn (sh' :. n) -			= rofu !: (sh' :. (k * n) `mod` len)+			= rofu ! (sh' :. (k * n) `mod` len)  	  in  A.sumAll $ A.zipWith (*) arrX wroots 

@@ -2,30 +2,32 @@ -- | 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-	( centerVector-	, centerMatrix)+	( center1d+	, center2d+	, center3d) where import Data.Array.Repa import Data.Array.Repa.Algorithms.Complex - -- | Apply the centering transform to a vector.-centerVector-	:: Array DIM1 Complex-	-> Array DIM1 Complex--{-# INLINE centerVector #-}-centerVector arr+center1d :: Array DIM1 Complex -> Array DIM1 Complex+{-# INLINE center1d #-}+center1d arr  = traverse arr id 	(\get ix@(_ :. x) -> ((-1) ^ x) * get ix)   -- | Apply the centering transform to a matrix.-centerMatrix-	:: Array DIM2 Complex-	-> Array DIM2 Complex--{-# INLINE centerMatrix #-}-centerMatrix arr+center2d :: Array DIM2 Complex -> Array DIM2 Complex+{-# INLINE center2d #-}+center2d arr  = traverse arr id 	(\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix)+++-- | Apply the centering transform to a 3d array.+center3d :: Array DIM3 Complex -> Array DIM3 Complex+{-# INLINE center3d #-}+center3d arr+ = traverse arr id+	(\get ix@(_ :. z :. y :. x) -> ((-1) ^ (z + y + x)) * get ix)

@@ -19,8 +19,9 @@  = force $ fromFunction sh f  where     f :: Shape sh => (sh :. Int) -> Complex-    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))-		:*: (- sin  (2 * pi * (fromIntegral i) / len))+    f (_ :. i) +	= ( cos  (2 * pi * (fromIntegral i) / len)+	  , - sin  (2 * pi * (fromIntegral i) / len))      len	= fromIntegral n @@ -36,7 +37,8 @@  = force $ fromFunction sh f  where     f :: Shape sh => (sh :. Int) -> Complex-    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))-		:*:   (sin  (2 * pi * (fromIntegral i) / len))+    f (_ :. i) +	= ( cos  (2 * pi * (fromIntegral i) / len)+	  , sin  (2 * pi * (fromIntegral i) / len))      len	= fromIntegral n

@@ -1,212 +1,185 @@-{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes #-}+{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes, ScopedTypeVariables, BangPatterns, FlexibleContexts #-}+{-# OPTIONS -fno-warn-incomplete-patterns #-} --- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm.------   Time complexity is O(n log n) in the size of the input.------   Input dimensions must be powers of two, else `error`.------   The `fft` and `ifft` functions (and friends) also compute the roots of unity needed.---   If you need to transform several arrays with the same extent then it is faster to---   compute the roots once using `calcRootsOfUnity` or `calcInverseRootsOfUnity`, ---   then call `fftWithRoots` directly.+-- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm. +--   Time complexity is O(n log n) in the size of the input.  -----   The inverse transforms provided also perform post-scaling so that `ifft` is the true inverse of `fft`. ---   If you don't want that then call `fftWithRoots` directly.+--   This uses a naive divide-and-conquer algorithm, the absolute performance is about+--   50x slower than FFTW in estimate mode. -----   The functions `fft2d` and `fft3d` require their inputs to be squares (and cubes) respectively. ---   This allows them to reuse the same roots-of-unity when transforming along each axis. If you ---   need to transform rectanglular arrays then call `fftWithRoots` directly. module Data.Array.Repa.Algorithms.FFT-	( fft,   ifft-	, fft2d, ifft2d-	, fft3d, ifft3d-	, fftWithRoots )+	( Mode(..)+	, isPowerOfTwo+	, fft3d+	, fft2d+	, fft1d) where-import Data.Array.Repa.Algorithms.DFT.Roots import Data.Array.Repa.Algorithms.Complex import Data.Array.Repa				as A-import Data.Ratio --- Vector Transform ---------------------------------------------------------------------------------- | Compute the DFT along the low order dimension of an array.-fft	:: Shape sh-	=> Array (sh :. Int) Complex-	-> Array (sh :. Int) Complex+data Mode+	= Forward+	| Reverse+	| Inverse+	deriving (Show, Eq) -fft v- = let	rofu	= calcRootsOfUnity (extent v)-   in	force $ fftWithRoots rofu v+{-# INLINE signOfMode #-}+signOfMode :: Mode -> Double+signOfMode mode+ = case mode of+	Forward		-> (-1)+	Reverse		->   1+	Inverse		->   1  --- | Compute the inverse DFT along the low order dimension of an array.-ifft	:: Shape sh-	=> Array (sh :. Int) Complex-	-> Array (sh :. Int) Complex+{-# INLINE isPowerOfTwo #-}+-- | Check if an `Int` is a power of two.+isPowerOfTwo :: Int -> Bool+isPowerOfTwo x+ = let	r	= (log (fromIntegral x) / log 2) :: Double+   in	ceiling r == (floor r :: Int) -ifft v- = let	_ :. len	= extent v-	scale		= fromIntegral len :*: 0-	rofu		= calcInverseRootsOfUnity (extent v)-   in	force $ A.map (/ scale) $ fftWithRoots rofu v  --- Matrix Transform ---------------------------------------------------------------------------------- | Compute the DFT of a square matrix.---   If the matrix is not square then `error`.-fft2d 	:: Array DIM2 Complex-	-> Array DIM2 Complex -fft2d arr- 	| Z :. height :. width	<- extent arr- 	, height /= width	-	= error $ "fft2d: height of matrix (" ++ show height ++ ")"-		++  " does not match width (" ++ show width  ++ ")"--	| otherwise-	= let	rofu		= calcRootsOfUnity (extent arr)-  		fftTrans 	= transpose . fftWithRoots rofu-   	  in	force $ fftTrans $ fftTrans arr----- | Compute the inverse DFT of a square matrix. -ifft2d	:: Array DIM2 Complex-	-> Array DIM2 Complex-	-ifft2d arr- 	| Z :. height :. width	<- extent arr- 	, height /= width	-	= error $ "fft2d: height of matrix (" ++ show height ++ ")"-		++  " does not match width (" ++ show width  ++ ")"+-- 3D Transform -----------------------------------------------------------------------------------+-- | Compute the DFT of a 3d array. Array dimensions must be powers of two else `error`.+fft3d 	:: Mode+	-> Array DIM3 Complex+	-> Array DIM3 Complex -	| otherwise-	= let	_ :. height :. width = extent arr-		scale		= fromIntegral (height * width) :*: 0-		rofu		= calcInverseRootsOfUnity (extent arr)-		fftTrans	= transpose . fftWithRoots rofu-	  in	force $ A.map (/ scale) $ fftTrans $ fftTrans arr-	+fft3d 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 "Data.Array.Repa.Algorithms.FFT: fft3d -- array dimensions must be powers of two."+	 else arr `deepSeqArray` +		case mode of+			Forward	-> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr+			Reverse	-> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr+			Inverse	-> force $ A.map (/ scale) +					$ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr --- Cube Transform ------------------------------------------------------------------------------------ | Compute the DFT of a 3d cube.---   If the array is not a cube then `error`.-fft3d 	:: Array DIM3 Complex+fftTrans3d +	:: Double+	-> Array DIM3 Complex  	-> Array DIM3 Complex -fft3d arrIn- 	| Z :. depth :. height :. width	<- extent arrIn- 	, (height /= width) || (height /= depth)-	= error $ "fft3d: array is not a cube"+{-# NOINLINE fftTrans3d #-}+fftTrans3d sign arr'+ = let 	arr		= force arr'+	(sh :. len)	= extent arr+   in	force $ rotate3d $ fft sign sh len arr -	| otherwise-	= let	rofu		= calcRootsOfUnity (extent arrIn) -		transpose3 arr-	 	 = traverse arr -        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) -            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) +rotate3d :: Array DIM3 Complex -> Array DIM3 Complex+{-# INLINE rotate3d #-}+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' -		fftTrans	= transpose3 . fftWithRoots rofu-	-  	  in	force $ fftTrans $ fftTrans $ fftTrans arrIn  --- | Compute the inverse DFT of a 3d cube.---   If the array is not a cube then `error`.-ifft3d 	:: Array DIM3 Complex-	-> Array DIM3 Complex+-- Matrix Transform -------------------------------------------------------------------------------+-- | Compute the DFT of a matrix. Array dimensions must be powers of two else `error`.+fft2d 	:: Mode+	-> Array DIM2 Complex+	-> Array DIM2 Complex -ifft3d arrIn- 	| Z :. depth :. height :. width	<- extent arrIn- 	, (height /= width) || (height /= depth)-	= error $ "ifft3d: array is not a cube"+fft2d mode arr+ = let	_ :. height :. width	= extent arr+	sign	= signOfMode mode+	scale 	= fromIntegral (width * height) +		+   in	if not (isPowerOfTwo height && isPowerOfTwo width)+	 then error "Data.Array.Repa.Algorithms.FFT: fft2d -- array dimensions must be powers of two."+	 else arr `deepSeqArray` +		case mode of+			Forward	-> fftTrans2d sign $ fftTrans2d sign arr+			Reverse	-> fftTrans2d sign $ fftTrans2d sign arr+			Inverse	-> force $ A.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr -	| otherwise-	= let	rofu		= calcInverseRootsOfUnity (extent arrIn)+fftTrans2d +	:: Double+	-> Array DIM2 Complex +	-> Array DIM2 Complex -		transpose3 arr-	 	 = traverse arr -        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) -            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) +{-# NOINLINE fftTrans2d #-}+fftTrans2d sign arr'+ = let 	arr		= force arr'+	(sh :. len)	= extent arr+   in	force $ transpose $ fft sign sh len arr -		_ :. depth :. height :. width -				= extent arrIn-		scale		= fromIntegral (height * width * depth) :*: 0 -		fftTrans	= transpose3 . fftWithRoots rofu-	  in	force $ A.map (/ scale) $ fftTrans $ fftTrans $ fftTrans arrIn--	--- Worker -------------------------------------------------------------------------------------------- | Generic function for computation of forward or inverse Discrete Fourier Transforms.---	Computation is along the low order dimension of the array.-fftWithRoots	-	:: forall sh-	.  Shape sh-	=> Array (sh :. Int) Complex		-- ^ Roots of unity.-	-> Array (sh :. Int) Complex		-- ^ Input values.-        -> Array (sh :. Int) Complex--fftWithRoots rofu v-	| not $ (denominator $ toRational (logBase (2 :: Double) $ fromIntegral vLen)) == 1-	= error $ "fft: vector length of " ++ show vLen ++ " is not a power of 2"+-- Vector Transform -------------------------------------------------------------------------------+-- | Compute the DFT of a vector. Array dimensions must be powers of two else `error`.+fft1d	:: Mode +	-> Array DIM1 Complex +	-> Array DIM1 Complex 	-	| rLen /= vLen-	= error $  "fft: length of vector (" ++ show vLen ++ ")"-		++ " does not match the length of the roots (" ++ show rLen ++ ")"+fft1d mode arr+ = let	_ :. len	= extent arr+	sign	= signOfMode mode+	scale	= fromIntegral len 	-	| otherwise-	= fftWithRoots' rofu v+   in	if not $ isPowerOfTwo len+	 then error "Data.Array.Repa.Algorithms.FFT: fft1d -- array dimensions must be powers of two."+	 else arr `deepSeqArray`+		case mode of+			Forward	-> fftTrans1d sign arr+			Reverse	-> fftTrans1d sign arr+			Inverse -> force $ A.map (/ scale) $ fftTrans1d sign arr -	where	_ :. rLen	= extent rofu-		_ :. vLen	= extent v+fftTrans1d+	:: Double +	-> Array DIM1 Complex+	-> Array DIM1 Complex -fftWithRoots'-	:: Shape sh-	=> Array (sh :. Int) Complex-	-> Array (sh :. Int) Complex-        -> Array (sh :. Int) Complex+{-# NOINLINE fftTrans1d #-}+fftTrans1d sign arr'+ = let	arr		= force arr'+	(sh :. len)	= extent arr+   in	fft sign sh len arr -{-# INLINE fftWithRoots' #-}-fftWithRoots' rofu v- = case extent v of-	_ :. 2	-> fft_two   v-	_	-> fft_split rofu v -{-# INLINE fft_two #-}-fft_two v- = let	vFn' vFn (sh :. 0)  = vFn (sh :. 0) + vFn (sh :. 1)-	vFn' vFn (sh :. 1)  = vFn (sh :. 0) - vFn (sh :. 1)-	vFn' _   _          = error "Data.Array.Repa.Algorithms.FFT fft_two fail"-   in	traverse v id vFn'+-- Rank Generalised Worker ------------------------------------------------------------------------+{-# INLINE fft #-}+fft !sign !sh !lenVec !vec+ = go lenVec 0 1+ where	go !len !offset !stride+	 | len == 2+	 = force $ fromFunction (sh :. 2) swivel 	-{-# INLINE fft_split #-}-fft_split rofu v- = let 	fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)+	 | otherwise+	 = combine len +		(go (len `div` 2) offset            (stride * 2))+		(go (len `div` 2) (offset + stride) (stride * 2)) -	fft_l  = traverse2 fft_lr rofu - 		   (\(sh :. 2 :. n) _ -> sh :. n)-		   (\f r (sh :. i)    -> f (sh :. 0 :. i) + r (sh :. i) * f (sh :. 1 :. i))+	 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))) -	fft_r  = traverse2 fft_lr rofu -		   (\(sh :. 2 :. n) _ -> sh :. n)-		   (\f r (sh :. i)    -> f (sh :. 0 :. i) - r (sh :. i) * f (sh :. 1 :. i))+		{-# INLINE combine #-}+		combine !len' 	evens@(Array _ [Region RangeAll GenManifest{}]) +				 odds@(Array _ [Region RangeAll GenManifest{}])+ 	 	 = evens `deepSeqArray` odds `deepSeqArray`+   	   	   let	odds'	= unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) +   	   	   in	force 	$ (evens +^ odds') A.++ (evens -^ odds') -   in	fft_l +:+ fft_r -{-# INLINE splitRofu #-}-splitRofu rofu- = traverse rofu-	(\(rSh :. rLen) 	-> rSh :. (2::Int) :. (rLen `div` 2))-	(\rFn (sh :. _ :. i) 	-> rFn (sh :. 2*i))--{-# INLINE splitVector #-}-splitVector v - = let	vFn' vFn (sh :. 0 :. i) = vFn (sh :. 2*i)-	vFn' vFn (sh :. 1 :. i) = vFn (sh :. 2*i+1)-	vFn' _   _              = error "Data.Array.Repa.Algorithms.FFT splitVector fail"+-- Compute a twiddle factor.+twiddle :: Double+	-> Int 			-- index+	-> Int 			-- length+	-> Complex -   in	traverse v-		(\(vSh :. vLen)    -> vSh :. 2 :. (vLen `div` 2)) -		vFn'-        +{-# INLINE twiddle #-}+twiddle sign k' n'+ 	=  (cos (2 * pi * k / n), sign * sin  (2 * pi * k / n))+	where 	k	= fromIntegral k'+		n	= fromIntegral n'+      

@@ -0,0 +1,45 @@+{-# LANGUAGE BangPatterns, PackageImports #-}+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}+module Data.Array.Repa.Algorithms.Iterate+	(iterateBlockwise, iterateBlockwise')+where+import Data.Array.Repa++-- | Iterate array transformation function a fixed number of times, applying `force2` between+--   each iteration.+iterateBlockwise+	:: (Elt a, Num a)+	=> Int					-- ^ Number of iterations to run for.+	-> (Array DIM2 a -> Array DIM2 a)	-- ^ Fn to step the array.+	-> Array DIM2 a				-- ^ Initial array value.+	-> Array DIM2 a+	+{-# INLINE iterateBlockwise #-}+iterateBlockwise steps f arrInit@(Array shInit [Region RangeAll (GenManifest vecInit)])+ = arrInit `deepSeqArray`+  goSolve steps shInit vecInit++ where	-- NOTE: We manually unpack the current array into its shape and vector to+	--	 stop GHC from unboxing the vector again for every loop. deepSeqing+	--	 the arrays at the start of solveLaplace makes the unboxings happen+	--	 at that point in the corresponding core code.+	goSolve !i !shCurrent !vecCurrent+	 = let	!arrCurrent	= fromVector shCurrent vecCurrent+	   in   if i == 0 +		 then arrCurrent+		 else let arrNew@(Array _ [Region RangeAll (GenManifest _)]) +				= force2 $ f arrCurrent+		      in  goSolve (i - 1) (extent arrNew) (toVector arrNew)+++-- | As above, but with the parameters flipped.+iterateBlockwise'+	:: (Elt a, Num a)+	=> Int					-- ^ Number of iterations to run for.+	-> Array DIM2 a				-- ^ Initial array value.+	-> (Array DIM2 a -> Array DIM2 a)	-- ^ Fn to step the array.+	-> Array DIM2 a++{-# INLINE iterateBlockwise' #-}+iterateBlockwise' steps arr fn+	= iterateBlockwise steps fn arr

@@ -1,4 +1,5 @@ {-# OPTIONS -fno-warn-incomplete-patterns #-}+{-# LANGUAGE PackageImports #-}  -- | Algorithms operating on matrices. -- @@ -11,20 +12,30 @@ module Data.Array.Repa.Algorithms.Matrix 	(multiplyMM) where-import Data.Array.Repa-	+import Data.Array.Repa	as A  -- | Matrix-matrix multiply.----multiplyMM	+multiplyMM 	:: Array DIM2 Double 	-> Array DIM2 Double 	-> Array DIM2 Double -multiplyMM  arr1 arr2- = multiplyMM' (force arr1) (force arr2)- where multiplyMM' arr1'@Manifest{} arr2'@Manifest{}-	= fold (+) 0 -  	$ traverse2 arr1' (force $ transpose arr2')-      		(\(sh :. m1 :. n1) -> \(_ :. n2 :. _m2) -> (sh :. m1 :. n2 :. n1))-		(\f1 -> \f2 -> \(sh :. i :. j :. k) -> f1 (sh :. i :. k) * f2 (sh :. j :. k))+{-# NOINLINE multiplyMM #-}+multiplyMM arr@(Array _ [Region RangeAll (GenManifest _)])+	   brr@(Array _ [Region RangeAll (GenManifest _)])+ = [arr, brr] `deepSeqArrays`+   A.force $ A.sum (A.zipWith (*) arrRepl brrRepl)+ where	trr@(Array _ [Region RangeAll (GenManifest _)])+			= force $ transpose2D brr+	arrRepl		= trr `deepSeqArray` A.extend (Z :. All   :. colsB :. All) arr+	brrRepl		= trr `deepSeqArray` A.extend (Z :. rowsA :. All   :. All) trr+	(Z :. _     :. rowsA) = extent arr+	(Z :. colsB :. _    ) = extent brr+	++transpose2D :: Elt e => Array DIM2 e -> Array DIM2 e+{-# INLINE transpose2D #-}+transpose2D arr+ = backpermute new_extent swap arr+ where	swap (Z :. i :. j)	= Z :. j :. i+	new_extent		= swap (extent arr)

@@ -0,0 +1,115 @@+{-# LANGUAGE BangPatterns #-}++module Data.Array.Repa.Algorithms.Randomish+ 	( 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+++-- | Use the ''minimal standard'' Lehmer generator to quickly generate some random+--   numbers with reasonable statistical properties. By ''reasonable'' we mean good+--   enough for games and test data, but not cryptography or anything where the+--   quality of the randomness really matters. +--+--   By nature of the algorithm, the maximum value in the output is clipped+--   to (valMin + 2^31 - 1)+-- +--   From ''Random Number Generators: Good ones are hard to find''+--   Stephen K. Park and Keith W. Miller.+--   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 sh Int		-- ^ Array of randomish numbers.++randomishIntArray !sh !valMin !valMax !seed+	= fromVector 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.++randomishIntVector !len !valMin' !valMax' !seed'+ = 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++	-- 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++	{-# INLINE f #-}+	f x		= multiplier * x `mod` modulus+ in G.create +     $ 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'++	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 sh Double	-- ^ Array of randomish numbers.++randomishDoubleArray !sh !valMin !valMax !seed+	= fromVector 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.++randomishDoubleVector !len !valMin !valMax !seed+ = let	range	= valMax - valMin++	mx	= 2^(30 :: Integer) - 1+	mxf	= fromIntegral mx+	ints	= randomishIntVector len 0 mx seed+	+   in	V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints

@@ -1,5 +1,5 @@ Name:                repa-algorithms-Version:             1.1.0.0+Version:             2.0.0.1 License:             BSD3 License-file:        LICENSE Author:              The DPH Team@@ -11,27 +11,33 @@ Homepage:            http://trac.haskell.org/repa Bug-reports:         http://trac.haskell.org/repa/newticket Description:-        NOTE: You must use the GHC head branch > 6.13.20100309 to get decent performance.         Reusable algorithms using the Repa array library.  Synopsis:         Algorithms using the Repa array library. -Tested-with: GHC == 6.13.20100309, GHC == 6.12.1- Library   Build-Depends:          base                 == 4.*,-        dph-base             == 0.4.*,-        repa                 == 1.1.*+        vector               == 0.7.*,+        repa                 == 2.0.*    ghc-options:-        -Odph -Wall -fno-warn-missing-signatures+        -Wall -fno-warn-missing-signatures+        -Odph+        -fsimplifier-phases=4+        -fstrictness-before=5+        -funfolding-use-threshold=30+        -funbox-strict-fields+        -fcpr-off    Exposed-modules:         Data.Array.Repa.Algorithms.Complex+        Data.Array.Repa.Algorithms.Randomish         Data.Array.Repa.Algorithms.DFT         Data.Array.Repa.Algorithms.DFT.Roots         Data.Array.Repa.Algorithms.DFT.Center         Data.Array.Repa.Algorithms.FFT         Data.Array.Repa.Algorithms.Matrix+        Data.Array.Repa.Algorithms.Convolve+        Data.Array.Repa.Algorithms.Iterate