repa-algorithms (empty) → 1.1.0.0

raw patch · 9 files changed

+531/−0 lines, 9 filesdep +basedep +dph-basedep +repasetup-changed

Dependencies added: base, dph-base, repa

Files

• Data/Array/Repa/Algorithms/Complex.hs +53/−0
• Data/Array/Repa/Algorithms/DFT.hs +100/−0
• Data/Array/Repa/Algorithms/DFT/Center.hs +31/−0
• Data/Array/Repa/Algorithms/DFT/Roots.hs +42/−0
• Data/Array/Repa/Algorithms/FFT.hs +212/−0
• Data/Array/Repa/Algorithms/Matrix.hs +30/−0
• LICENSE +24/−0
• Setup.hs +2/−0
• repa-algorithms.cabal +37/−0

@@ -0,0 +1,53 @@+{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}++-- | Strict complex doubles.+module Data.Array.Repa.Algorithms.Complex+	( Complex+	, mag+	, arg+	, (:*:)(..))+where+import 	Data.Array.Parallel.Base ((:*:)(..))++-- | Strict complex doubles.+type Complex +	= 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+++instance Fractional Complex where+  (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+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)+ = normaliseAngle $ atan2 im re++ where 	normaliseAngle :: Double -> Double+	normaliseAngle f+	 | f < - pi	+	 = normaliseAngle (f + 2 * pi)+	+	 | f > pi+	 = normaliseAngle (f - 2 * pi)++	 | otherwise+	 = f

@@ -0,0 +1,100 @@+{-# LANGUAGE TypeOperators, RankNTypes, PatternGuards #-}++-- | Compute the Discrete Fourier Transform (DFT) along the low order dimension+--   of an array. +--+--   This uses the naive algorithm and takes O(n^2) time. +--   However, you can transform an array with an arbitray extent, unlike with FFT which requires+--   each dimension to be a power of two.+--+--   The `dft` and `idft` functions 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 `dftWithRoots` directly.+--+--   You can also compute single values of the transform using `dftWithRootsSingle`.+module Data.Array.Repa.Algorithms.DFT +	( dft+	, idft+	, dftWithRoots+	, dftWithRootsSingle)+where+import Data.Array.Repa.Algorithms.DFT.Roots+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa				as A++-- | Compute the DFT along the low order dimension of an array.+dft 	:: forall sh+	.  Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++dft v+ = let	rofu	= calcRootsOfUnity (extent v)+   in	force $ dftWithRoots rofu v+++-- | Compute the inverse DFT along the low order dimension of an array.+idft 	:: forall sh+	.  Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++idft v+ = let	_ :. len	= extent v+	scale		= fromIntegral len :*: 0+	rofu		= calcInverseRootsOfUnity (extent v)+   in	force $ A.map (/ scale) $ dftWithRoots rofu v+++-- | 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`.+dftWithRoots+	:: forall sh+	.  Shape sh+	=> Array (sh :. Int) Complex		-- ^ Roots of unity.+	-> Array (sh :. Int) Complex		-- ^ Input array.+	-> Array (sh :. Int) Complex++dftWithRoots rofu arr+	| _ :. 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 ++ ")"++	| otherwise+	= traverse arr id (\_ k -> dftWithRootsSingle rofu arr k)+		++-- | Compute a single value of the DFT.+--	The extent of the given roots must match that of the input array, else `error`.+dftWithRootsSingle+	:: forall sh+	.  Shape sh+	=> Array (sh :. Int) Complex 		-- ^ Roots of unity.+	-> Array (sh :. Int) Complex		-- ^ Input array.+	-> (sh :. Int)				-- ^ Index of the value we want.+	-> Complex++{-# INLINE dftWithRootsSingle #-}+dftWithRootsSingle rofu arrX (_ :. k)+	| _ :. 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 ++ ")"++	| 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)++	  in  A.sumAll $ A.zipWith (*) arrX wroots++

@@ -0,0 +1,31 @@++-- | 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)+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+ = 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+ = traverse arr id+	(\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix)

@@ -0,0 +1,42 @@+{-# LANGUAGE TypeOperators, RankNTypes #-}++-- | Calculation of roots of unity for the forward and inverse DFT\/FFT.+module Data.Array.Repa.Algorithms.DFT.Roots+	( calcRootsOfUnity+	, calcInverseRootsOfUnity)+where+import Data.Array.Repa+import Data.Array.Repa.Algorithms.Complex++-- | Calculate roots of unity for the forward transform.+calcRootsOfUnity+	:: forall sh+	.  Shape sh+	=> (sh :. Int) 			-- ^ Length of lowest dimension of result.+	-> Array (sh :. Int) Complex++calcRootsOfUnity sh@(_ :. n) + = 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))++    len	= fromIntegral n+++-- | Calculate roots of unity for the inverse transform.+calcInverseRootsOfUnity+	:: forall sh+	.  Shape sh+	=> (sh :. Int) 			-- ^ Length of lowest dimension of result.+	-> Array (sh :. Int) Complex++calcInverseRootsOfUnity sh@(_ :. n) + = 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))++    len	= fromIntegral n

@@ -0,0 +1,212 @@+{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes #-}++-- | 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.+--+--   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.+--+--   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 )+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++fft v+ = let	rofu	= calcRootsOfUnity (extent v)+   in	force $ fftWithRoots rofu v+++-- | Compute the inverse DFT along the low order dimension of an array.+ifft	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++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  ++ ")"++	| 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+	++-- Cube Transform ---------------------------------------------------------------------------------+-- | Compute the DFT of a 3d cube.+--   If the array is not a cube then `error`.+fft3d 	:: 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"++	| 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)) ++		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++ifft3d arrIn+ 	| Z :. depth :. height :. width	<- extent arrIn+ 	, (height /= width) || (height /= depth)+	= error $ "ifft3d: array is not a cube"++	| otherwise+	= let	rofu		= calcInverseRootsOfUnity (extent arrIn)++		transpose3 arr+	 	 = traverse arr +        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) +            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) ++		_ :. 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"+	+	| rLen /= vLen+	= error $  "fft: length of vector (" ++ show vLen ++ ")"+		++ " does not match the length of the roots (" ++ show rLen ++ ")"+	+	| otherwise+	= fftWithRoots' rofu v++	where	_ :. rLen	= extent rofu+		_ :. vLen	= extent v++fftWithRoots'+	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex+        -> Array (sh :. Int) Complex++{-# 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'+	+{-# INLINE fft_split #-}+fft_split rofu v+ = let 	fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)++	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))++	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))++   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"++   in	traverse v+		(\(vSh :. vLen)    -> vSh :. 2 :. (vLen `div` 2)) +		vFn'+        

@@ -0,0 +1,30 @@+{-# OPTIONS -fno-warn-incomplete-patterns #-}++-- | Algorithms operating on matrices.+-- +--   These functions should give performance comparable with nested loop C+--   implementations, but not block-based, cache friendly, SIMD using, vendor+--   optimised implementions. +--   If you care deeply about runtime performance then you may be better off using +--   a binding to LAPACK, such as hvector.+--+module Data.Array.Repa.Algorithms.Matrix+	(multiplyMM)+where+import Data.Array.Repa+	++-- | Matrix-matrix multiply.+--+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))

@@ -0,0 +1,24 @@+Copyright (c) 2010, University of New South Wales.+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:+    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.+    * Redistributions in binary form must reproduce the above copyright+      notice, this list of conditions and the following disclaimer in the+      documentation and/or other materials provided with the distribution.+    * Neither the name of the University of New South Wales nor the+      names of its contributors may be used to endorse or promote products+      derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,37 @@+Name:                repa-algorithms+Version:             1.1.0.0+License:             BSD3+License-file:        LICENSE+Author:              The DPH Team+Maintainer:          Ben Lippmeier <benl@ouroborus.net>+Build-Type:          Simple+Cabal-Version:       >=1.6+Stability:           experimental+Category:            Data Structures+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.*++  ghc-options:+        -Odph -Wall -fno-warn-missing-signatures++  Exposed-modules:+        Data.Array.Repa.Algorithms.Complex+        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