packages feed

vector-fftw 0.1.3.8 → 0.1.4.0

raw patch · 9 files changed

+480/−17 lines, 9 filesdep +QuickCheckdep +test-frameworkdep +test-framework-quickcheck2dep ~basedep ~primitivedep ~vector

Dependencies added: QuickCheck, test-framework, test-framework-quickcheck2, vector-fftw

Dependency ranges changed: base, primitive, vector

Files

+ Changelog.md view
@@ -0,0 +1,12 @@+# vector-fftw changelog++## 0.1.4.0++* Introduce multi-dimensional transforms:+  * `Numeric.FFT.Vector.Unitary.Multi`+  * `Numeric.FFT.Vector.Invertible.Multi`+  * `Numeric.FFT.Vector.Unnormalized.Multi`++## 0.1 through 0.1.3.8++The pre-historic era.
Numeric/FFT/Vector/Base.hsc view
@@ -6,6 +6,11 @@             PlanType(..),             plan,             run,+            -- * multi-demensional Transforms+            TransformND(..),+            planOfTypeND,+            planND,+            runND,             -- * Plans             Plan(),             planInputSize,@@ -36,7 +41,7 @@ import Control.Monad(forM_) import Foreign (Storable(..), Ptr, FunPtr,                 ForeignPtr, withForeignPtr, newForeignPtr)-import Foreign.C (CInt(..), CUInt)+import Foreign.C (CInt(..), CUInt, CSize(..)) import Data.Bits ( (.|.) ) import Data.Complex(Complex(..)) import Foreign.Storable.Complex()@@ -154,7 +159,7 @@ ------------------ -- Malloc/free of fftw array -foreign import ccall unsafe fftw_malloc :: CInt -> IO (Ptr a)+foreign import ccall unsafe fftw_malloc :: CSize -> IO (Ptr a) foreign import ccall "&" fftw_free :: FunPtr (Ptr a -> IO ())  newFFTVector :: forall a . Storable a => Int -> IO (MS.MVector RealWorld a)@@ -186,14 +191,14 @@     planInput <- newFFTVector m_in     planOutput <- newFFTVector m_out     MS.unsafeWith planInput $ \inP -> MS.unsafeWith planOutput $ \outP -> do-    pPlan <- makePlan (toEnum n) inP outP $ planInitFlags ptype DestroyInput-    cPlan <- newPlan pPlan-    -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs-    -- aren't released too soon:-    let planExecute = MS.unsafeWith planInput $ \_ ->-                        MS.unsafeWith planOutput $ \_ ->-                          withPlan cPlan fftw_execute-    return $ normalization n $ Plan {..}+        pPlan <- makePlan (toEnum n) inP outP $ planInitFlags ptype DestroyInput+        cPlan <- newPlan pPlan+        -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs+        -- aren't released too soon:+        let planExecute = MS.unsafeWith planInput $ \_ ->+                            MS.unsafeWith planOutput $ \_ ->+                                withPlan cPlan fftw_execute+        return $ normalization n $ Plan {..}   where     m_in = inputSize n     m_out = outputSize n@@ -212,6 +217,67 @@             (planOfType Estimate p $ creationSizeFromInput p $ V.length v)             v {-# INLINE run #-}++-----------------------+-- TransformND: methods of plan creation for multi-dimensional plans.++-- | A transform which may be applied to vectors of different sizes.+--+-- @since 0.2+data TransformND a b = TransformND {+                        inputSizeND :: Int -> Int,+                        outputSizeND :: Int -> Int,+                        creationSizeFromInputND :: Int -> Int,+                        makePlanND :: CInt -> Ptr CInt -> Ptr a -> Ptr b -> CFlags -> IO (Ptr CPlan),+                        normalizationND :: VS.Vector Int -> Plan a b -> Plan a b+                    }++-- | Create a 'Plan' of a specific size for this transform.+-- 'dims' must have rank greater or equal to 1+--+-- @since 0.2+planOfTypeND :: (Storable a, Storable b) => PlanType+                                -> TransformND a b -> VS.Vector Int -> Plan a b+planOfTypeND ptype TransformND{..} dims+  | m_in <= 0 || m_out <= 0 = error "Can't (yet) plan for empty arrays!"+  | otherwise  = unsafePerformIO $ do+    mdims <- unsafeThaw $ V.map toEnum dims+    planInput <- newFFTVector m_in+    planOutput <- newFFTVector m_out+    MS.unsafeWith mdims $ \dimsP -> MS.unsafeWith planInput $ \inP -> MS.unsafeWith planOutput $ \outP -> do+        pPlan <- makePlanND (toEnum $ V.length dims) dimsP inP outP $ planInitFlags ptype DestroyInput+        cPlan <- newPlan pPlan+        -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs+        -- aren't released too soon:+        let planExecute = MS.unsafeWith planInput $ \_ ->+                                MS.unsafeWith planOutput $ \_ ->+                                withPlan cPlan fftw_execute+        return $ normalizationND dims $ Plan {..}+  where+    m = V.product $ V.init dims+    m_in = m * inputSizeND (V.last dims)+    m_out = m * outputSizeND (V.last dims)+{-# INLINE planOfTypeND #-}++-- | Create a 'Plan' of a specific size.  This function is equivalent to+-- @'planOfType' 'Estimate'@.+--+-- @since 0.2+planND :: (Storable a, Storable b) => TransformND a b -> VS.Vector Int -> Plan a b+planND = planOfTypeND Estimate+{-# INLINE planND #-}++-- | Create and run a 'Plan' for the given transform.+--+-- @since 0.2+runND :: (Vector v a, Vector v b, Storable a, Storable b)+            => TransformND a b -> VS.Vector Int ->  v a -> v b+runND p = \dims v ->+  let creationSize = V.init dims `V.snoc` creationSizeFromInputND p (V.last dims) in+    execute+    (planOfTypeND Estimate p creationSize)+    v+{-# INLINE runND #-}  --------------------------- -- For scaling input/output:
+ Numeric/FFT/Vector/Invertible/Multi.hs view
@@ -0,0 +1,47 @@+{- |+This module provides  normalized multi-dimensional versions of the transforms in @fftw@.++The forwards transforms in this module are identical to those in "Numeric.FFT.Vector.Unnormalized".+The backwards transforms are normalized to be their inverse operations (approximately, due to floating point precision).++For more information on the underlying transforms, see+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.++@since 0.2+-}++module Numeric.FFT.Vector.Invertible.Multi+  (+        -- * Creating and executing 'Plan's+        run,+        plan,+        execute,+        -- * Complex-to-complex transforms+        U.dft,+        idft,+        -- * Real-to-complex transforms+        U.dftR2C,+        dftC2R,+  ) where++import Numeric.FFT.Vector.Base+import qualified Numeric.FFT.Vector.Unnormalized.Multi as U+import Data.Complex+import qualified Data.Vector.Storable as VS++-- | A backward discrete Fourier transform which is the inverse of 'U.dft'.  The output and input sizes are the same (@n@).+idft :: TransformND (Complex Double) (Complex Double)+idft = U.idft {normalizationND = \ns -> constMultOutput $ 1 / toEnum (VS.product ns)}++-- | A normalized backward discrete Fourier transform which is the left inverse of+-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)+--+-- This 'Transform' behaves differently than the others:+--+--  - Calling @planND dftC2R dims@ where @dims = [n0, ..., nk]@ creates a 'Plan' whose /output/ size is @dims@, and whose+--    /input/ size is @[n0, ..., nk \`div\` 2 + 1]@.+--+--  - If @length v == n0 * ... * nk@, then @length (run dftC2R v) == n0 * ... * 2*(nk-1)@.+--+dftC2R :: TransformND (Complex Double) Double+dftC2R = U.dftC2R {normalizationND = \ns -> constMultOutput $ 1 / toEnum (VS.product ns)}
Numeric/FFT/Vector/Plan.hs view
@@ -5,6 +5,10 @@                 PlanType(..),                 plan,                 run,+                TransformND(),+                planOfTypeND,+                planND,+                runND,                 -- * Plans                 Plan(),                 planInputSize,
+ Numeric/FFT/Vector/Unitary/Multi.hs view
@@ -0,0 +1,89 @@+{- |+This module provides normalized versions of the transforms in @fftw@.++All of the transforms are normalized so that++ - Each transform is unitary, i.e., preserves the inner product and the sum-of-squares norm of its input.++ - Each backwards transform is the inverse of the corresponding forwards transform.++(Both conditions only hold approximately, due to floating point precision.)++For more information on the underlying transforms, see+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.+--+-- @since 0.2+-}++module Numeric.FFT.Vector.Unitary.Multi+  (+        -- * Creating and executing 'Plan's+        run,+        plan,+        execute,+        -- * Complex-to-complex transforms+        dft,+        idft,+        -- * Real-to-complex transforms+        dftR2C,+        dftC2R,+  ) where++import Control.Exception (assert)+import Control.Monad (forM_)+import Numeric.FFT.Vector.Base+import qualified Numeric.FFT.Vector.Unnormalized.Multi as U+import Data.Complex+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as MS+import Control.Monad.Primitive(RealWorld)++-- | A discrete Fourier transform. The output and input sizes are the same (@n@).+--+-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(-2pi i j k\/n)@+dft :: TransformND (Complex Double) (Complex Double)+dft = U.dft {normalizationND = \ns -> constMultOutput $ 1 / sqrt (toEnum (VS.product ns))}++-- | An inverse discrete Fourier transform.  The output and input sizes are the same (@n@).+--+-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(2pi i j k\/n)@+idft :: TransformND (Complex Double) (Complex Double)+idft = U.idft {normalizationND = \ns -> constMultOutput $ 1 / sqrt (toEnum (VS.product ns))}++-- | A forward discrete Fourier transform with real data.  If the input size is @n@,+-- the output size will be @n \`div\` 2 + 1@.+dftR2C :: TransformND Double (Complex Double)+dftR2C = U.dftR2C {normalizationND = \ns -> modifyOutput $+                    complexR2CScaling (sqrt 2) ns (outputSizeND U.dftR2C $ VS.last ns)+        }++-- | A normalized backward discrete Fourier transform which is the left inverse of+-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)+--+-- This 'Transform' behaves differently than the others:+--+--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose+--    /input/ size is @n \`div\` 2 + 1@.+--+--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.+--+dftC2R :: TransformND (Complex Double) Double+dftC2R = U.dftC2R {normalizationND = \ns -> modifyInput $+                    complexR2CScaling (sqrt 0.5) ns (inputSizeND U.dftC2R $ VS.last ns)+        }++complexR2CScaling :: Double -> VS.Vector Int -> Int -> MS.MVector RealWorld (Complex Double) -> IO ()+complexR2CScaling !t !ns !len !a = assert (MS.length a == VS.product (VS.init ns) * len) $ do+    let !s1 = sqrt (1/toEnum (VS.product ns))+    let !s2 = t * s1+    -- Justification for the use of unsafeModify:+    -- The output size is 2n+1; so if n>0 then the output size is >=1;+    -- and if n even then the output size is >=3.+    forM_ [0.. VS.product (VS.init ns) - 1] $ \idx -> do+      unsafeModify a (idx * len) $ scaleByD s1+      if odd (VS.last ns)+        then multC s2 (MS.unsafeSlice (idx * len + 1) (len-1) a)+        else do+            unsafeModify a (idx * len + len - 1) $ scaleByD s1+            multC s2 (MS.unsafeSlice (idx * len + 1) (len-2) a)+
Numeric/FFT/Vector/Unnormalized.hsc view
@@ -107,6 +107,7 @@             normalization = const id         } + r2rTransform :: CKind -> Transform Double Double r2rTransform kind = Transform {                     inputSize = id,
+ Numeric/FFT/Vector/Unnormalized/Multi.hsc view
@@ -0,0 +1,92 @@+{- |+Raw, unnormalized multi-dimensional versions of the transforms in @fftw@.++Note that the forwards and backwards transforms of this module are not actually+inverses.  For example, @run idft (run dft v) /= v@ in general.++For more information on the individual transforms, see+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.++@since 0.2+-}++module Numeric.FFT.Vector.Unnormalized.Multi+  (+        -- * Creating and executing 'Plan's+        run,+        plan,+        execute,+        -- * Complex-to-complex transforms+        dft,+        idft,+        -- * Real-to-complex transforms+        dftR2C,+        dftC2R,+  ) where++import Numeric.FFT.Vector.Base+import Foreign+import Foreign.C+import Data.Complex++#include <fftw3.h>++-- | Whether the complex fft is forwards or backwards.+type CDirection = CInt++foreign import ccall unsafe fftw_plan_dft+    :: CInt -> Ptr CInt -> Ptr (Complex Double) -> Ptr (Complex Double)+        -> CDirection -> CFlags -> IO (Ptr CPlan)++foreign import ccall unsafe fftw_plan_dft_r2c+    :: CInt -> Ptr CInt -> Ptr Double -> Ptr (Complex Double) -> CFlags+        -> IO (Ptr CPlan)++foreign import ccall unsafe fftw_plan_dft_c2r+    :: CInt -> Ptr CInt -> Ptr (Complex Double) -> Ptr Double -> CFlags+        -> IO (Ptr CPlan)++dftND :: CDirection -> TransformND (Complex Double) (Complex Double)+dftND d = TransformND+  { inputSizeND = id+  , outputSizeND = id+  , creationSizeFromInputND = id+  , makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft rk dims a b d+  , normalizationND = const id+  }++-- | A forward discrete Fourier transform.  The output and input sizes are the same (@n@).+dft :: TransformND (Complex Double) (Complex Double)+dft = dftND (#const FFTW_FORWARD)++-- | A backward discrete Fourier transform.  The output and input sizes are the same (@n@).+idft :: TransformND (Complex Double) (Complex Double)+idft = dftND (#const FFTW_BACKWARD)++-- | A forward discrete Fourier transform with real data.  If the input size is @n0 * ... * nk@,+-- the output size will be @n0 * ... * nk \`div\` 2 + 1@.+dftR2C :: TransformND Double (Complex Double)+dftR2C = TransformND {+              inputSizeND = id,+              outputSizeND = \n -> n `div` 2 + 1,+              creationSizeFromInputND = id,+              makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft_r2c rk dims a b,+              normalizationND = const id+          }++-- | A backward discrete Fourier transform which produces real data.+--+-- This 'Transform' behaves differently than the others:+--+--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose+--    /input/ size is @n \`div\` 2 + 1@.+--+--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.+dftC2R :: TransformND (Complex Double) Double+dftC2R = TransformND {+            inputSizeND = \n -> n `div` 2 + 1,+            outputSizeND = id,+            creationSizeFromInputND = \n -> 2 * (n-1),+            makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft_c2r rk dims a b,+            normalizationND = const id+        }
+ tests/FFTProperties.hs view
@@ -0,0 +1,130 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+-- This module uses the test-framework-quickcheck2 package.+module Main where++import Control.Monad+import qualified Data.Vector.Unboxed as V+import qualified Data.Vector.Storable as VS+import Data.Complex++import Test.Framework (defaultMain, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.QuickCheck++import qualified Numeric.FFT.Vector.Invertible as I+import qualified Numeric.FFT.Vector.Invertible.Multi as IM+import qualified Numeric.FFT.Vector.Unitary as U+import qualified Numeric.FFT.Vector.Unitary.Multi as UM+import Numeric.FFT.Vector.Plan++main = defaultMain+            -- NB: There's no explicit tests for the Unnormalized package.+            -- However, its Planners are implicitly used by the other modules,+            -- so it's covered in the below tests.+            [ testGroup "invertibility"+              [ testProperty "I.dft" $ prop_invert I.dft I.idft+              , testProperty "I.dftR2C" $ prop_invert I.dftR2C I.dftC2R+              , testProperty "I.dct1" $ prop_invert I.dct1 I.idct1+              , testProperty "I.dct2" $ prop_invert I.dct2 I.idct2+              , testProperty "I.dct3" $ prop_invert I.dct3 I.idct3+              , testProperty "I.dct4" $ prop_invert I.dct4 I.idct4+              , testProperty "I.dst1" $ prop_invert I.dst1 I.idst1+              , testProperty "I.dst2" $ prop_invert I.dst2 I.idst2+              , testProperty "I.dst3" $ prop_invert I.dst3 I.idst3+              , testProperty "I.dst4" $ prop_invert I.dst4 I.idst4+              , testProperty "U.dft" $ prop_invert U.dft U.idft+              , testProperty "U.dftR2C" $ prop_invert U.dftR2C U.dftC2R+              , testProperty "U.dct2" $ prop_invert U.dct2 U.idct2+              ]+            , testGroup "orthogonality"+              [ testProperty "U.dft" $ prop_orthog U.dft+              , testProperty "U.idft" $ prop_orthog U.idft+              , testProperty "U.dftR2C" $ prop_orthog U.dftR2C+              , testProperty "U.dftC2R" $ prop_orthog U.dftR2C+              , testProperty "U.dct2" $ prop_orthog U.dct2+              , testProperty "U.idct2" $ prop_orthog U.idct2+              , testProperty "U.dct4" $ prop_orthog U.dct4+              ]+            , testGroup "invertibility ND"+              [ testProperty "IM.dft" $ prop_invertND IM.dft IM.idft+              , testProperty "IM.dftR2C" $ prop_invertND IM.dftR2C IM.dftC2R+              , testProperty "UM.dft" $ prop_invertND UM.dft UM.idft+              , testProperty "UM.dftR2C" $ prop_invertND UM.dftR2C UM.dftC2R+              ]+            , testGroup "orthogonality"+              [ testProperty "UM.dft" $ prop_orthogND UM.dft+              , testProperty "UM.idft" $ prop_orthogND UM.idft+              , testProperty "UM.dftR2C" $ prop_orthogND UM.dftR2C+              , testProperty "UM.dftC2R" $ prop_orthogND UM.dftR2C+              ]+            ]++-------------------+-- An instance of Arbitrary that probably belongs in another package.++instance (V.Unbox a, Arbitrary a) => Arbitrary (V.Vector a) where+    arbitrary = V.fromList `fmap` arbitrary+++-------------------------+-- Support functions to compare Doubles for (near) equality.++class Num a => Mag a where+    mag :: a -> Double++instance Mag Double where+    mag = abs++instance Mag (Complex Double) where+    mag = magnitude++-- Robustly test whether two Doubles are nearly identical.+close :: Mag a => a -> a -> Bool+close x y = tol > mag (x-y) / max 1 (mag x + mag y)+  where+    tol = 1e-10++withinTol :: (Mag a, V.Unbox a) => V.Vector a -> V.Vector a -> Bool+withinTol a b+    | V.length a /= V.length b = False+    | otherwise = V.and $ V.zipWith close a b+++---------------------+-- The actual properties++-- Test whether the inverse actually inverts the forward transform.+prop_invert f g a = let+                        p1 = plan f (V.length a)+                        p2 = plan g (V.length a)+                    in (V.length a > 1) ==> withinTol a $ execute p2 $ execute p1 a++-- Test whether the transform preserves the L2 (sum-of-squares) norm.+prop_orthog f a = let+                    p1 = plan f (V.length a)+                  in (V.length a > 1) ==> close (norm2 a) (norm2 $ execute p1 a)++data DimsAndValues a = DimsAndValues (VS.Vector Int) (V.Vector a)+  deriving (Show)++instance (Arbitrary a, V.Unbox a) => Arbitrary (DimsAndValues a) where+  arbitrary = do+    dims <- liftM (VS.fromList . map getPositive) arbitrary `suchThatMap` maybeReduceSize+    values <- V.replicateM (VS.product dims) arbitrary+    return (DimsAndValues dims values)+    where+      -- We use this to prevent test cases from growing too big+      maybeReduceSize ds =+        if VS.product ds < 1000 then Just ds else maybeReduceSize (VS.init ds)++prop_invertND f g (DimsAndValues ds a) = let+                        p1 = planND f ds+                        p2 = planND g ds+                    in (V.length a > 1) ==> withinTol a $ execute p2 $ execute p1 a++prop_orthogND f (DimsAndValues ds a) = let+                    p1 = planND f ds+                  in (V.length a > 1) ==> close (norm2 a) (norm2 $ execute p1 a)++norm2 a = sqrt $ V.sum $ V.map (\x -> x*x) $ V.map mag a
vector-fftw.cabal view
@@ -1,6 +1,8 @@+cabal-version:       >=1.10+ Name:                vector-fftw -Version:             0.1.3.8+Version:             0.1.4.0 License:             BSD3 License-file:        LICENSE Author:              Judah Jacobson@@ -8,13 +10,12 @@ Copyright:           (c) Judah Jacobson, 2010 Category:            Math Build-type:          Simple-Cabal-version:       >=1.6 Homepage:            http://hackage.haskell.org/package/vector-fftw Synopsis:            A binding to the fftw library for one-dimensional vectors. Description:         This package provides bindings to the fftw library for one-dimensional vectors.                      It provides both high-level functions and more low-level manipulation of fftw plans.                      .-                     We provide three different modules which wrap fftw's operations:+                     We provide three different modules which wrap @fftw@'s operations:                      .                       - "Numeric.FFT.Vector.Unnormalized" contains the raw transforms;                      .@@ -22,7 +23,11 @@                      .                       - "Numeric.FFT.Vector.Unitary" additionally scales all transforms to preserve the L2 (sum-of-squares) norm of the                         input.-Tested-With:         GHC == 7.6.2, GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.1+                     .+                     In addition, we provide @.Multi@ modules for each of these providing multi-dimensional+                     transforms.+Extra-Source-Files:  Changelog.md+Tested-With:         GHC == 7.6.2, GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.1, GHC == 8.6.4, GHC == 8.8.3, GHC == 8.10.1  source-repository head     type:   git@@ -30,26 +35,43 @@   Library+  default-language: Haskell2010   Exposed-modules:         Numeric.FFT.Vector.Unnormalized+        Numeric.FFT.Vector.Unnormalized.Multi         Numeric.FFT.Vector.Invertible+        Numeric.FFT.Vector.Invertible.Multi         Numeric.FFT.Vector.Unitary+        Numeric.FFT.Vector.Unitary.Multi         Numeric.FFT.Vector.Plan    Other-modules:         Numeric.FFT.Vector.Base -  Build-depends: base>=4.3 && < 4.11,+  Build-depends: base>=4.3 && < 4.15,                  vector>=0.9 && < 0.13,-                 primitive>=0.6 && < 0.7,+                 primitive>=0.6 && < 0.8,                  storable-complex==0.2.*   if os(windows)     Extra-libraries: fftw3-3   else     Extra-libraries: fftw3 -  Extensions: ForeignFunctionInterface, RecordWildCards, BangPatterns, FlexibleInstances,+  default-extensions: ForeignFunctionInterface, RecordWildCards, BangPatterns, FlexibleInstances,                 ScopedTypeVariables   ghc-options: -Wall    Ghc-Options: -O2++test-suite properties+  default-language: Haskell2010+  ghc-options: -Wall -threaded+  type: exitcode-stdio-1.0+  hs-source-dirs: tests+  main-is: FFTProperties.hs+  build-depends: base,+                 QuickCheck,+                 test-framework,+                 test-framework-quickcheck2,+                 vector,+                 vector-fftw