repa-fftw (empty) → 3.2.3.1
raw patch · 6 files changed
+427/−0 lines, 6 filesdep +basedep +carraydep +criterionsetup-changed
Dependencies added: base, carray, criterion, deepseq, fft, random, repa, repa-algorithms, repa-fftw, storable-complex, tasty, tasty-hunit, tasty-quickcheck
Files
- Data/Array/Repa/FFTW.hs +180/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- exec/bench.hs +49/−0
- exec/test.hs +94/−0
- repa-fftw.cabal +72/−0
+ Data/Array/Repa/FFTW.hs view
@@ -0,0 +1,180 @@+{-|+Module : $Header$+CopyRight : (c) 8c6794b6, 2011, 2012+License : BSD3+Maintainer : 8c6794b6@gmail.com+Stability : experimental+Portability : portable++Performs fft of repa array data via FFTW.++Currently supporting ('Complex' Double) arrays for dimensions 'DIM1', 'DIM2',+and 'DIM3' only.++-}++module Data.Array.Repa.FFTW+ ( -- * Examples+ -- $examples++ -- * Multi dimension functions+ -- $multi++ -- * References+ -- $references++ -- * FFT functions (1 dimension)+ fft+ , ifft+ -- * FFT functions (2 dimension)+ , fft2d+ , ifft2d+ -- * FFT functions (3 dimension)+ , fft3d+ , ifft3d+ ) where++import Data.Complex (Complex(..))+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Storable (Storable(..))+import System.IO.Unsafe (unsafePerformIO)++import Data.Array.CArray (CArray)+import Data.Array.Repa ((:.)(..), Array, DIM1, DIM2, DIM3, Z(..))+import Data.Array.Repa.Repr.ForeignPtr (F)+import Foreign.Storable.Complex ()++import qualified Data.Array.CArray as C+import qualified Data.Array.Repa as R+import qualified Data.Array.Repa.Repr.ForeignPtr as RF+import qualified Math.FFT as FFT++{-$examples++Sample module:++> import Data.Complex+> import Data.Array.Repa+> import Data.Array.Repa.Eval+> import Data.Array.Repa.Repr.ForeignPtr+> import Data.Array.Repa.FFTW+>+> a :: Array F DIM1 (Complex Double)+> a = fromList (Z :. 4) [i :+ 0 | i <- [0..3]]++Loading above in ghci:++>>> toList a+[0.0 :+ 0.0,1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]+>>> toList $ fft a+[6.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]+>>> toList $ ifft $ fft a+[0.0 :+ 0.0,1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]++-}++{-$multi++Although FFTW library has more flexiblity, choice of dimensions in+multi-dimension FFT functions ('fft2d', 'ifft2d' ... etc) were using fixed+dimensions to perform FFTs. Those choices were made after @fft2dP@, @fft3dP@+functions from @repa-algorithms@ package.++-}++{-$references++* fftw : <http://www.fftw.org>++* fftw haskell binding : <http://hackage.haskell.org/package/fft>++* repa-algorithms: <http://hackage.haskell.org/package/repa-algorithms>++-}++-- --------------------------------------------------------------------------+-- Exposed functions++-- | Performs 1 dimension forward fft.+fft :: Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double)+fft = c2r . FFT.dft . r2c+{-# INLINE fft #-}++-- | Performs 1 dimension inverse fft.+ifft :: Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double)+ifft = c2r . FFT.idft . r2c+{-# INLINE ifft #-}++-- | Performs 2 dimension forward fft.+fft2d :: Array F DIM2 (Complex Double) -> Array F DIM2 (Complex Double)+fft2d = c2r2d . FFT.dftN [0,1] . r2c2d+{-# INLINE fft2d #-}++-- | Performs 2 dimension inverse fft.+ifft2d :: Array F DIM2 (Complex Double) -> Array F DIM2 (Complex Double)+ifft2d = c2r2d . FFT.idftN [0,1] . r2c2d+{-# INLINE ifft2d #-}++-- | Performs 3 dimension forward fft.+fft3d :: Array F DIM3 (Complex Double) -> Array F DIM3 (Complex Double)+fft3d = c2r3d . FFT.dftN [0,1,2] . r2c3d+{-# INLINE fft3d #-}++-- | Performs 3 dimension inverse fft.+ifft3d :: Array F DIM3 (Complex Double) -> Array F DIM3 (Complex Double)+ifft3d = c2r3d . FFT.idftN [0,1,2] . r2c3d+{-# INLINE ifft3d #-}++-- --------------------------------------------------------------------------+-- Guts++r2c :: Array F DIM1 (Complex Double) -> CArray Int (Complex Double)+r2c rarr = unsafePerformIO $ do+ let _:.nelem = R.extent rarr+ fptr = RF.toForeignPtr rarr+ C.unsafeForeignPtrToCArray fptr (0,nelem-1)+{-# INLINE r2c #-}++c2r :: CArray Int (Complex Double) -> Array F DIM1 (Complex Double)+c2r carr = case C.toForeignPtr carr of+ (n, fptr) -> let sh = Z:.n in+ R.computeS $ R.fromFunction sh $ \ix ->+ unsafePerformIO $ withForeignPtr fptr $ \ptr ->+ peekElemOff ptr $ R.toIndex sh ix+{-# INLINE c2r #-}++r2c2d :: Array F DIM2 (Complex Double) -> CArray (Int, Int) (Complex Double)+r2c2d rarr = unsafePerformIO $ do+ let _:.n1:.n2 = R.extent rarr+ fptr = RF.toForeignPtr rarr+ C.unsafeForeignPtrToCArray fptr ((0,0), (n1-1, n2-1))+{-# INLINE r2c2d #-}++c2r2d :: CArray (Int, Int) (Complex Double) -> Array F DIM2 (Complex Double)+c2r2d carr = case C.toForeignPtr carr of+ (n, fptr) ->+ let sh = Z:.n':.n'+ n' = ceiling $ (sqrt $ fromIntegral n :: Double)+ in R.computeS $ R.fromFunction sh $ \ix ->+ unsafePerformIO $ withForeignPtr fptr $ \ptr ->+ peekElemOff ptr $ R.toIndex sh ix+{-# INLINE c2r2d #-}++r2c3d :: Array F DIM3 (Complex Double)+ -> CArray (Int, Int, Int) (Complex Double)+r2c3d rarr = unsafePerformIO $ do+ let _:.n1:.n2:.n3 = R.extent rarr+ fptr = RF.toForeignPtr rarr+ C.unsafeForeignPtrToCArray fptr ((0,0,0), (n1-1, n2-1, n3-1))+{-# INLINE r2c3d #-}++c2r3d :: CArray (Int, Int, Int) (Complex Double)+ -> Array F DIM3 (Complex Double)+c2r3d carr = case C.toForeignPtr carr of+ (n, fptr) ->+ let sh = Z:.n':.n':.n'+ n' = ceiling $ fromIntegral n ** (1/3 :: Double)+ in R.computeS $ R.fromFunction sh $ \ix ->+ unsafePerformIO $ withForeignPtr fptr $ \ptr ->+ peekElemOff ptr $ R.toIndex sh ix+{-# INLINE c2r3d #-}
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2011-2013, 8c6794b6++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 8c6794b6 nor the names of other+ 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 AND CONTRIBUTORS+"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 THE COPYRIGHT+OWNER OR CONTRIBUTORS 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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ exec/bench.hs view
@@ -0,0 +1,49 @@+{-# LANGUAGE FlexibleInstances #-}+{-|+Module : $Header$+CopyRight : (c) 8c6794b6, 2011, 2012+License : BSD3+Maintainer : 8c6794b6@gmail.com+Stability : experimental+Portability : portable++Benchmark for comparing FFT with repa-fftw to repa-algorithms.+-}+module Main where++import Control.DeepSeq (NFData(..))+import Data.Complex+import System.Random++import Criterion.Main++import Data.Array.Repa ((:.)(..), Array, U, DIM1, Z(..))+import Data.Array.Repa.Repr.ForeignPtr (F)+import qualified Data.Array.Repa as R+import qualified Data.Array.Repa.Eval as Eval+import qualified Data.Array.Repa.Algorithms.FFT as FFT+import qualified Data.Array.Repa.FFTW as FFTW++type RFFTArr = Array U DIM1 (Double, Double)++main :: IO ()+main = do+ rs <- randomRs (0,1) `fmap` newStdGen+ is <- randomRs (0,1) `fmap` newStdGen+ let bench_fft n =+ let mkarr ks = Eval.fromList (Z:.n) $ take n ks+ ts = R.fromListUnboxed (Z:.n) $ take n $ zip rs is+ cs = mkarr $ zipWith (:+) rs is+ ra = FFT.fft1dP FFT.Forward :: RFFTArr -> IO RFFTArr+ in ts `seq` cs `seq`+ bgroup ("n="++show n)+ [ bench "repa-algorithms" (nfIO (ra ts))+ , bench "repa-fftw" (nf FFTW.fft cs)+ ]+ defaultMain $ map (\k -> bench_fft (2^k)) [9..13::Int]++instance NFData (Array U DIM1 (Double, Double)) where+ rnf arr = arr `R.deepSeqArray` ()++instance NFData (Array F DIM1 (Complex Double)) where+ rnf arr = arr `R.deepSeqArray` ()
+ exec/test.hs view
@@ -0,0 +1,94 @@+-- | Main entry point for test executable.++module Main where++import Data.Array.Repa ((:.)(..), Array, DIM1, DIM2, DIM3, Shape, Z(..), toList)+import Data.Array.Repa.Eval (fromList)+import Data.Array.Repa.Repr.ForeignPtr (F)+import Data.Complex (Complex(..), magnitude)++import Test.Tasty (defaultMain, testGroup)+import Test.Tasty.HUnit (Assertion, assertBool, testCase)+import Test.Tasty.QuickCheck (Property, choose, forAll, testProperty, vector)++import Data.Array.Repa.FFTW (fft, ifft, fft2d, ifft2d, fft3d, ifft3d)+++main :: IO ()+main = defaultMain $+ testGroup "repa-fftw"+ [ testCase "magsum1d-hunit" case_magsum1d+ , testCase "magsum2d-hunit" case_magsum2d+ , testCase "magsum3d-hunit" case_magsum3d+ , testProperty "magsum1d-qc" prop_magsum1d+ , testProperty "magsum2d-qc" prop_magsum2d+ , testProperty "magsum3d-qc" prop_magsum3d+ ]++case_magsum1d :: Assertion+case_magsum1d = do+ let a = fromList (Z :. 4) [i :+ 0 | i <- [0..3]]+ assertBool "diff is greater than epsilon" $ compareMagnitudeSum1d a++case_magsum2d :: Assertion+case_magsum2d = do+ let a = fromList (Z :. 4 :. 4) [i :+ 0 | i <- [0..15]]+ assertBool "diff is greater than epsilon" $ compareMagnitudeSum2d a++case_magsum3d :: Assertion+case_magsum3d = do+ let a = fromList (Z :. 4 :. 4 :. 4) [i :+ 0 | i <- [0..63]]+ assertBool "diff is greater than epsilon" $ compareMagnitudeSum3d a++prop_magsum1d :: Property+prop_magsum1d =+ forAll (choose (1, 16::Int)) $ \pow ->+ let size = (2 ^ pow)+ in forAll (vector size) $ \vals ->+ compareMagnitudeSum1d $ fromList (Z :. size) vals++prop_magsum2d :: Property+prop_magsum2d =+ forAll (choose (1, 4::Int)) $ \pow ->+ let extnt = 2 ^ pow+ size = extnt ^ (2 :: Int)+ in forAll (vector size) $ \vals ->+ compareMagnitudeSum2d $ fromList (Z :. extnt :. extnt) vals++prop_magsum3d :: Property+prop_magsum3d =+ forAll (choose (1, 4::Int)) $ \pow ->+ let extnt = 2 ^ pow+ size = extnt ^ (3 :: Int)+ in forAll (vector size) $ \vals ->+ compareMagnitudeSum3d $ fromList (Z :. extnt :. extnt :. extnt) vals++-- | Performs fft and ifft, then compareing magnitude of sum with the+-- original.+compareMagnitudeSum1d :: Array F DIM1 (Complex Double) -> Bool+compareMagnitudeSum1d = compareMagnitudeSum fft ifft++-- | Like 'compareMagnitudeSum', for 2-dimensional array.+compareMagnitudeSum2d :: Array F DIM2 (Complex Double) -> Bool+compareMagnitudeSum2d = compareMagnitudeSum fft2d ifft2d++-- | Like 'compareMagnitudeSum', for 3-dimensional array.+compareMagnitudeSum3d :: Array F DIM3 (Complex Double) -> Bool+compareMagnitudeSum3d = compareMagnitudeSum fft3d ifft3d++-- | Apply given functions, then compare the result with original.+compareMagnitudeSum ::+ (Shape sh)+ => (Array F sh (Complex Double) -> a)+ -- ^ Forward FFT function.+ -> (a -> Array F sh (Complex Double))+ -- ^ Inverse FFT function.+ -> Array F sh (Complex Double)+ -- ^ Original input.+ -> Bool+compareMagnitudeSum fn ifn arr = diff < epsilon+ where+ diff = abs (mag arr - mag arr')+ mag = magnitude . sum . toList+ arr' = ifn $ fn arr+ epsilon = 1e-2
+ repa-fftw.cabal view
@@ -0,0 +1,72 @@+Name: repa-fftw+Version: 3.2.3.1+Synopsis: Perform fft with repa via FFTW+License: BSD3+License-file: LICENSE+Author: <8c6794b6@gmail.com>+Maintainer: <8c6794b6@gmail.com>+Category: Math+Build-type: Simple+Cabal-version: >=1.8++Description:+ Performs FFT on repa array data with fftw.++flag dev+ description:+ Activate developer flags+ default:+ False++source-repository head+ type: git+ location: https://github.com/8c6794b6/repa-fftw.git++Library+ ghc-options:+ -Wall -fno-warn-orphans+ -O3 -fllvm -optl-O3+ exposed-modules:+ Data.Array.Repa.FFTW+ build-depends:+ base >= 4.6 && < 5,+ carray >= 0.1.5 && < 0.2,+ fft >= 0.1.5 && < 0.2,+ repa >= 3.2.3,+ storable-complex >= 0.2++test-suite test+ type:+ exitcode-stdio-1.0+ hs-source-dirs:+ exec+ main-is:+ test.hs+ ghc-options:+ -Wall+ build-depends:+ base >= 4.6 && < 5,+ repa >= 3.2.3,+ repa-fftw -any,+ tasty >= 0.4.0,+ tasty-hunit >= 0.4.1,+ tasty-quickcheck >= 0.3.1++benchmark bench+ type:+ exitcode-stdio-1.0+ hs-source-dirs:+ exec+ main-is:+ bench.hs+ ghc-options:+ -Wall -fno-warn-orphans+ -threaded -rtsopts+ build-depends:+ base >= 4.6 && < 5,+ criterion >= 0.8,+ deepseq >= 1.3.0 && < 2,+ repa >= 3.2.3,+ repa-algorithms >= 3.2,+ repa-fftw -any,+ random >= 1.0.1.1