{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} module Test.Data.Array.Accelerate.Fourier where -- (tests) where import qualified Data.Array.Accelerate.Fourier.Real as FourierReal import qualified Data.Array.Accelerate.Fourier.Preprocessed as Prep import qualified Data.Array.Accelerate.Fourier.Adhoc as Adhoc import qualified Data.Array.Accelerate.Fourier.Planned as Planned import qualified Data.Array.Accelerate.Convolution.Preprocessed as ConvPrep import qualified Data.Array.Accelerate.Convolution.Adhoc as Convolution import qualified Data.Array.Accelerate.Cyclic as Cyclic import qualified Data.Array.Accelerate.Interpreter as AI import qualified Data.Array.Accelerate as A import Data.Array.Accelerate.Fourier.Planned (Transform, ) import Data.Array.Accelerate.Fourier.Utility (scaleDown, ) import Data.Array.Accelerate (Acc, Exp, Array, DIM1, DIM2, DIM3, Z(Z), (:.)((:.)), ) import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg import qualified Data.Array.Accelerate.Utility.Sliced as Sliced import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp import Data.Array.Accelerate.Utility.Lift.Exp (expr) import qualified Test.QuickCheck as QC import Test.QuickCheck (Arbitrary, arbitrary, quickCheck, ) import Control.Monad (liftM2, liftM3, guard, ) import qualified Data.Array.Accelerate.Data.Complex as Complex import Data.Complex (Complex((:+)), cis, ) tolerance :: Double tolerance = 1e-10 approxEqualAbs :: (A.RealFloat a) => Exp a -> Exp a -> Exp a -> Exp Bool approxEqualAbs eps x y = abs (x-y) A.<= eps approxEqualComplexAbs :: (A.RealFloat a) => Exp a -> Exp (Complex a) -> Exp (Complex a) -> Exp Bool approxEqualComplexAbs eps x y = Complex.magnitude (x-y) A.<= eps genComplex :: QC.Gen (Complex Double) genComplex = liftM2 (:+) (QC.choose (-1,1)) (QC.choose (-1,1)) newtype Normed0 = Normed0 (Exp (Complex Double)) deriving (Show) instance Arbitrary Normed0 where arbitrary = fmap (Normed0 . A.constant) genComplex data Normed1 = Normed1 Int (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1 where arbitrary = fmap (\xs -> let len = length xs in Normed1 len $ A.use $ A.fromList (Z :. len) xs) $ QC.listOf genComplex floorPowerOfTwo :: Int -> Int floorPowerOfTwo len = 2 ^ (floor (logBase 2 (fromIntegral len :: Double)) :: Int) data Normed1PowerTwo = Normed1PowerTwo Int (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1PowerTwo where arbitrary = fmap (\xs -> let len = floorPowerOfTwo $ length xs in Normed1PowerTwo len $ A.use $ A.fromList (Z :. len) xs) $ liftM2 (:) genComplex (QC.listOf genComplex) {- For every reasonable pair of powers of 3 and 5 it computes the largest power of 2, such that their product is below @n@. -} floor5Smooth :: Int -> Int floor5Smooth n = fromInteger $ maximum $ map (maximum . floorSmooths 2 5 (fromIntegral n)) $ floorSmooths 2 3 (fromIntegral n) $ fromIntegral $ floorPowerOfTwo n {- | @floorSmooths a b n m@ replaces successively @a@ factors in @m@ by @b@ factors while keeping the product below @n@. -} floorSmooths :: Integer -> Integer -> Integer -> Integer -> [Integer] floorSmooths a b n = let divMany k = case divMod k a of (q,r) -> guard (r==0) >> if q>n then divMany q else go q go m = m : divMany (m*b) in go data Normed1Smooth5 = Normed1Smooth5 Int (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1Smooth5 where arbitrary = fmap (\xs -> let len = floor5Smooth $ length xs in Normed1Smooth5 len $ A.use $ A.fromList (Z :. len) xs) $ liftM2 (:) genComplex (QC.listOf genComplex) data Normed1Even = Normed1Even Int (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1Even where arbitrary = fmap (\xs -> let len = 2 * length xs in Normed1Even len $ A.use $ A.fromList (Z :. len) $ concatMap (\(x0,x1) -> [x0,x1]) xs) $ QC.listOf $ liftM2 (,) genComplex genComplex data Normed1Pair = Normed1Pair Int (Acc (Array DIM1 (Complex Double))) (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1Pair where arbitrary = fmap (\xys -> let len = length xys (xs, ys) = unzip xys in Normed1Pair len (A.use $ A.fromList (Z :. len) xs) (A.use $ A.fromList (Z :. len) ys)) $ QC.listOf $ liftM2 (,) genComplex genComplex data Normed1Triple = Normed1Triple Int (Acc (Array DIM1 (Complex Double))) (Acc (Array DIM1 (Complex Double))) (Acc (Array DIM1 (Complex Double))) deriving (Show) instance Arbitrary Normed1Triple where arbitrary = fmap (\xyzs -> let len = length xyzs (xs, ys, zs) = unzip3 xyzs in Normed1Triple len (A.use $ A.fromList (Z :. len) xs) (A.use $ A.fromList (Z :. len) ys) (A.use $ A.fromList (Z :. len) zs)) $ QC.listOf $ liftM3 (,,) genComplex genComplex genComplex data Normed2 = Normed2 Int Int (Acc (Array DIM2 (Complex Double))) deriving (Show) instance Arbitrary Normed2 where arbitrary = do xs <- liftM2 (:) QC.arbitrary $ QC.listOf genComplex let len = length xs height <- QC.choose (1, round (sqrt (fromIntegral len :: Double))) let width = div len height return $ Normed2 width height $ A.use $ A.fromList (Z :. height :. width) xs data Normed3 = Normed3 Int Int Int (Acc (Array DIM3 (Complex Double))) deriving (Show) instance Arbitrary Normed3 where arbitrary = do xs <- liftM2 (:) QC.arbitrary $ QC.listOf genComplex let len = length xs lenThd = round $ (fromIntegral len :: Double) ** recip 3 height <- QC.choose (1, lenThd) width <- QC.choose (1, lenThd) let depth = div len (width*height) return $ Normed3 width height depth $ A.use $ A.fromList (Z :. depth :. height :. width) xs -- duplicate of Private.cycleDim3 cycleDim3 :: (A.Elt a) => Transform DIM3 a cycleDim3 arr = A.backpermute (Exp.modify (expr:.expr:.expr:.expr) (\(sh:.k:.m:.n) -> (sh:.n:.k:.m)) $ A.shape arr) (Exp.modify (expr:.expr:.expr:.expr) (\(ix:.n:.k:.m) -> (ix:.k:.m:.n))) arr basisVector :: Int -> Int -> Array DIM1 (Complex Double) basisVector len freq = A.fromList (Z:.len) $ map (\k -> cis $ 2*pi * fromIntegral (k*freq) / fromIntegral len) $ iterate (1+) 0 norm2 :: Acc (Array DIM1 (Complex Double)) -> Acc (A.Scalar Double) norm2 = A.map sqrt . A.sum . A.map (Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i) scalarProduct :: Acc (Array DIM1 (Complex Double)) -> Acc (Array DIM1 (Complex Double)) -> Acc (A.Scalar (Complex Double)) scalarProduct xs ys = A.foldAll (+) 0 $ A.zipWith (*) xs (A.map Complex.conjugate ys) complexFromReal :: Acc (Array DIM1 Double) -> Acc (Array DIM1 (Complex Double)) complexFromReal = A.map (A.lift . (:+0)) toSelfAdjoint :: Transform DIM1 (Complex Double) toSelfAdjoint x = A.zipWith (+) x $ A.map Complex.conjugate $ Cyclic.reverse x infixl 6 .+. (.+.) :: Acc (Array DIM1 (Complex Double)) -> Acc (Array DIM1 (Complex Double)) -> Acc (Array DIM1 (Complex Double)) (.+.) = A.zipWith (+) quickCheckWithSign :: (Arbitrary normed, Show normed) => (Planned.Sign Double -> normed -> Bool) -> IO () quickCheckWithSign = quickCheck {- should be replaced by (A.==) in future -} class (A.Shape sh, A.Slice sh) => EqShape sh where eqShape :: Exp sh -> Exp sh -> Exp Bool instance EqShape Z where eqShape _ _ = A.constant True instance (EqShape sh, i ~ Int) => EqShape (sh:.i) where eqShape = Exp.modify2 (expr:.expr) (expr:.expr) $ \(sh0:.n0) (sh1:.n1) -> n0 A.== n1 A.&& eqShape sh0 sh1 infix 4 =~= (=~=) :: (EqShape sh) => Acc (Array sh (Complex Double)) -> Acc (Array sh (Complex Double)) -> Acc (A.Scalar Bool) (=~=) xs ys = A.map (eqShape (A.shape xs) (A.shape ys) A.&&) $ A.and $ A.flatten $ A.zipWith (approxEqualComplexAbs (A.constant tolerance)) xs ys run :: Acc (A.Scalar Bool) -> Bool run = Acc.the . AI.run tests :: [(String, IO ())] tests = ("fourier generic vs. preprocessed dit2", quickCheck $ \sign (Normed1PowerTwo len x) -> run $ Planned.transform sign len x =~= Prep.dit2 sign len x) : ("fourier generic vs. preprocessed dif2", quickCheck $ \sign (Normed1PowerTwo len x) -> run $ Planned.transform sign len x =~= Prep.dif2 sign len x) : ("fourier generic vs. preprocessed ditSplitRadix", quickCheck $ \sign (Normed1PowerTwo len x) -> run $ Planned.transform sign len x =~= Prep.ditSplitRadix sign len x) : ("fourier generic vs. adhoc dit2", quickCheck $ \sign (Normed1PowerTwo len x) -> run $ Planned.transform sign len x =~= Adhoc.dit2 (A.constant sign) x) : ("fourier generic vs. adhoc ditSplitRadix", quickCheck $ \sign (Normed1PowerTwo len x) -> run $ Planned.transform sign len x =~= Adhoc.ditSplitRadix (A.constant sign) x) : ("fourier generic vs. adhoc dit235", quickCheck $ \sign (Normed1Smooth5 len x) -> run $ Planned.transform sign len x =~= Adhoc.dit235 (A.constant sign) x) : ("fourier adhoc chirp 2 vs. chirp 235", quickCheck $ \sign (Normed1 _len x) -> run $ Adhoc.transformChirp2 (A.constant sign) x =~= Adhoc.transformChirp235 (A.constant sign) x) : ("fourier generic vs. adhoc auto", quickCheck $ \sign (Normed1 len x) -> run $ Planned.transform sign len x =~= Adhoc.transform (A.constant sign) x) : ("fourier generic vs. adhoc chirp 235", quickCheck $ \sign (Normed1 len x) -> run $ Planned.transform sign len x =~= Adhoc.transformChirp235 (A.constant sign) x) : ("fourier generic vs. chirp2", quickCheck $ \sign (Normed1 len x) -> run $ Planned.transform sign len x =~= Planned.transformChirp2 sign len x) : ("fourier generic vs. chirp235", quickCheck $ \sign (Normed1 len x) -> run $ Planned.transform sign len x =~= Planned.transformChirp235 sign len x) : ("homogeneity", quickCheck $ \sign (Normed0 x) (Normed1 len y) -> run $ let transform = Planned.transform sign len in transform (A.map (x*) y) =~= A.map (x*) (transform y)) : ("additivity", quickCheck $ \sign (Normed1Pair len x y) -> run $ let transform = Planned.transform sign len in A.zipWith (+) (transform x) (transform y) =~= transform (A.zipWith (+) x y)) : ("basis vector", quickCheck $ \(Normed1 len _x) kp -> run $ let transform = Planned.transform Planned.inverse len k = mod kp len unit = A.use $ A.fromList (Z:.len) $ replicate k 0 ++ 1 : repeat 0 in transform unit =~= A.use (basisVector len k)) : ("fourier inverse", quickCheck $ \(Normed1 len x) -> run $ x =~= (scaleDown $ Planned.transform Planned.forward len $ Planned.transform Planned.inverse len x)) : ("double fourier = reverse", quickCheck $ \sign (Normed1 len x) -> run $ let transform = Planned.transform sign len in x =~= (Cyclic.reverse $ scaleDown $ transform $ transform x)) : ("fourier of reverse", quickCheck $ \sign (Normed1 len x) -> run $ let transform = Planned.transform sign len in Cyclic.reverse (transform x) =~= transform (Cyclic.reverse x)) : ("fourier of conjugate", quickCheck $ \sign (Normed1 len x) -> run $ let transform = Planned.transform sign len in (A.map Complex.conjugate $ transform x) =~= (transform $ A.map Complex.conjugate $ Cyclic.reverse x)) : ("isometry", quickCheck $ \sign (Normed1 len x) -> run $ let transform = Planned.transform sign len in A.zipWith (approxEqualAbs $ A.constant tolerance) (norm2 $ transform x) (A.map (A.constant (sqrt (fromIntegral len)) *) $ norm2 x)) : ("unitarity", quickCheck $ \sign (Normed1Pair len x y) -> run $ let transform = Planned.transform sign len in A.zipWith (approxEqualComplexAbs $ A.constant tolerance) (scalarProduct (transform x) (transform y)) (A.map (A.constant (fromIntegral len) *) $ scalarProduct x y)) : ("convolution commutativity", quickCheck $ \(Normed1Pair len x y) -> run $ let (.*.) = Planned.convolveCyclic len in x .*. y =~= y .*. x) : ("convolution associativity", quickCheck $ \(Normed1Triple len x y z) -> run $ let (.*.) = Planned.convolveCyclic len in (x .*. y) .*. z =~= x .*. (y .*. z)) : ("convolution distributivity", quickCheck $ \(Normed1Triple len x y z) -> run $ let (.*.) = Planned.convolveCyclic len in x .*. (y .+. z) =~= (x .*. y) .+. (x .*. z)) : ("convolution karatsuba rec vs. loop", quickCheck $ \(Normed1 len xy) -> run $ let x = A.map Complex.real xy y = A.map Complex.imag xy in complexFromReal (ConvPrep.karatsuba len x y) =~= complexFromReal (Convolution.karatsuba x y)) : {- No instance for (A.IsNum (Complex Double)) arising from a use of 'ConvPrep.karatsuba' -} ("convolution karatsuba", quickCheck $ \(Normed1Pair len x y) -> run $ let resultLen = max 0 $ 2*len-1 in Convolution.complex Convolution.karatsuba x y =~= Planned.convolveCyclic resultLen (Sliced.pad 0 (A.constant resultLen) x) (Sliced.pad 0 (A.constant resultLen) y)) : ("convolution cyclic karatsuba", quickCheck $ \(Normed1Pair len x y) -> run $ Convolution.complex (Convolution.cyclic Convolution.karatsuba) x y =~= Planned.convolveCyclic len (Sliced.pad 0 (A.constant len) x) (Sliced.pad 0 (A.constant len) y)) : ("real to spectrum", quickCheck $ \(Normed1Even len x) -> run $ let xr = A.map Complex.real x in FourierReal.toSpectrum (Planned.transform Planned.forward (div len 2)) xr =~= Planned.transform Planned.forward len (complexFromReal xr)) : ("real from spectrum", quickCheck $ \(Normed1Even len x) -> run $ let xSelfAdjoint = toSelfAdjoint x in (complexFromReal $ FourierReal.fromSpectrum (Planned.transform Planned.inverse (div len 2)) xSelfAdjoint) =~= Planned.transform Planned.inverse len xSelfAdjoint) : ("real to and from spectrum", quickCheck $ \(Normed1Even len x) -> run $ let xr = A.map Complex.real x len2 = div len 2 in (scaleDown $ complexFromReal $ FourierReal.fromSpectrum (Planned.transform Planned.inverse len2) $ FourierReal.toSpectrum (Planned.transform Planned.forward len2) xr) =~= complexFromReal xr) : ("real from and to spectrum", quickCheck $ \(Normed1Even len x) -> run $ let len2 = div len 2 xSelfAdjoint = toSelfAdjoint x in (scaleDown $ FourierReal.toSpectrum (Planned.transform Planned.forward len2) $ FourierReal.fromSpectrum (Planned.transform Planned.inverse len2) $ xSelfAdjoint) =~= xSelfAdjoint) : ("double real to spectrum, even", quickCheck $ \(Normed1Even len x) -> run $ let xr = A.map Complex.real x xi = A.map Complex.imag x transform = Planned.transform Planned.forward (div len 2) (specr,speci) = A.unzip $ FourierReal.untangleSpectra $ Planned.transform Planned.forward len x in A.zipWith (A.&&) (FourierReal.toSpectrum transform xr =~= specr) (FourierReal.toSpectrum transform xi =~= speci)) : ("double real to spectrum, arbitrary", quickCheck $ \(Normed1 len x) -> run $ let xr = complexFromReal $ A.map Complex.real x xi = complexFromReal $ A.map Complex.imag x transform = Planned.transform Planned.forward len (specr,speci) = A.unzip $ FourierReal.untangleSpectra $ transform x in A.zipWith (A.&&) (transform xr =~= specr) (transform xi =~= speci)) : ("entangle and untangle spectrum of real data", quickCheck $ \(Normed1Pair _len x y) -> run $ let (xt,yt) = A.unzip $ A.map (A.uncurry FourierReal.untangleCoefficient) $ A.zipWith FourierReal.entangleCoefficient x y in A.zipWith (A.&&) (x =~= xt) (y =~= yt)) : ("double real from spectrum", quickCheck $ \(Normed1 len x) -> run $ let imagUnit = A.constant (0:+1) xSelfAdjoint = toSelfAdjoint x ySelfAdjoint = toSelfAdjoint $ A.map (imagUnit*) x transform = Planned.transform Planned.inverse len (xSignal,ySignal) = A.unzip $ FourierReal.twoFromSpectrum transform $ A.zip xSelfAdjoint ySelfAdjoint in A.zipWith (A.&&) (transform xSelfAdjoint =~= A.map (A.lift . (:+0)) xSignal) (transform ySelfAdjoint =~= A.map (A.lift . (:+0)) ySignal)) : ("transform2d vs. transposition, preprocessed", quickCheckWithSign $ \sign (Normed2 width height x) -> run $ let transformH = Prep.transform2d (Prep.SubTransformPair (Planned.transform sign width) (Planned.transform sign height)) transformV = Prep.transform2d (Prep.SubTransformPair (Planned.transform sign height) (Planned.transform sign width)) in LinAlg.transpose (transformH x) =~= transformV (LinAlg.transpose x)) : ("transform2d vs. transposition, adhoc", quickCheckWithSign $ \sign (Normed2 _width _height x) -> run $ let transform = Adhoc.transform2d (Adhoc.SubTransform (Adhoc.transformChirp2 (A.constant sign))) in LinAlg.transpose (transform x) =~= transform (LinAlg.transpose x)) : ("transform3d vs. transposition, preprocessed", quickCheckWithSign $ \sign (Normed3 width height depth x) -> run $ let transformH = Prep.transform3d (Prep.SubTransformTriple (Planned.transform sign width) (Planned.transform sign height) (Planned.transform sign depth)) transformV = Prep.transform3d (Prep.SubTransformTriple (Planned.transform sign height) (Planned.transform sign depth) (Planned.transform sign width)) in cycleDim3 (transformH x) =~= transformV (cycleDim3 x)) : ("transform2d vs. transposition, adhoc", quickCheckWithSign $ \sign (Normed3 _width _height _depth x) -> run $ let transform = Adhoc.transform2d (Adhoc.SubTransform (Adhoc.transformChirp2 (A.constant sign))) in LinAlg.transpose (transform x) =~= transform (LinAlg.transpose x)) : []