content: -- Test of the computational correctness of my current FFT candidate.
  --
  -- Original author: David Banas <capn.freako@gmail.com>
  -- Original date:   October 3, 2015
  --
  -- Copyright (c) 2015 David Banas; all rights reserved World wide.
  --
  -- NOTE: I have moved general FFT expression development up to the
  --       top level of the TreeViz project, into the fft-ccc.hs file.
  --
  --       Please, conduct all future explorations in general FFT
  --       expression, using Conal's machinery, which aren't necessarily
  --       intended to be compilable, but rather as a less constrained
  --       exploration of FFT expression alternatives, there.
  --
  --       I'm making this change of venue, for two reasons:
  --
  --       1) I find this file very useful, specifically, as a staging
  --          area for new test cases, intended for insertion into Conal's
  --          TreeTest.hs file, where they will actually be compiled
  --          by his machinery.
  --
  --       2) I'm realizing that a less constrained (by the requirements
  --          for compilation) vehicle for general exploration of more
  --          abstract, higher level expressions of the FFT operation,
  --          using Conal's data types and there associated functions, as
  --          well as others (potentially) would be very useful.
  --
  --       3) I'd like to start a wiki page, which would track the history
  --          of this more generic development of FFT expression, and the
  --          TreeViz project seems a natural base for such an effort,
  --          desparately needing a wiki of its own to narrate the history
  --          of the effort to search for more efficient FFT implementations.
  
  {-# LANGUAGE GADTs #-}
  {-# LANGUAGE TemplateHaskell #-}
  {-# LANGUAGE TypeSynonymInstances #-}
  {-# LANGUAGE MultiParamTypeClasses #-}
  {-# LANGUAGE FlexibleInstances #-}
  {-# LANGUAGE FlexibleContexts #-}
  
  module Main where
  
  import Prelude hiding ({- id,(.), -}foldl,foldr,sum,product,zipWith,reverse,and,or,scanl,minimum,maximum)
  
  import Control.Applicative
  import Control.Arrow
  import Control.Monad (forM_, unless)
  import Data.Complex
  import Data.Foldable (Foldable, sum, foldl', toList)
  import Data.Newtypes.PrettyDouble
  import System.Exit (exitFailure)
  import TypeUnary.Nat (IsNat, Nat(..), nat, N2, N3, N4, N5)  -- , N6)
  
  -- import Test.QuickCheck (choose, vectorOf, elements, collect)
  import Test.QuickCheck (choose, vectorOf)
  import Test.QuickCheck.Arbitrary
  import Test.QuickCheck.All (quickCheckAll)
  
  import Circat.Scan (lproducts, LScan)
  import qualified Circat.Pair as P
  import qualified Circat.RTree as RT
  import Circat.RTree (bottomSplit)
  
  type RTree = RT.Tree
  
  -- FFT, as a class
  -- (The LScan constraint comes from the use of 'lproducts', in 'addPhase'.)
  class (LScan f) => FFT f a where
      fft  :: f a -> f a  -- Computes the FFT of a functor.
  
  -- Note that this definition of the FFT instance for Pair assumes DIT.
  -- How can we eliminate this assumption and make this more general?
  instance (RealFloat a, Applicative f, Foldable f, Num (f (Complex a)), FFT f (Complex a)) => FFT P.Pair (f (Complex a)) where
      fft = P.inP (uncurry (+) &&& uncurry (-)) . P.secondP addPhase . fmap fft
  
  instance (IsNat n, RealFloat a) => FFT (RTree n) (Complex a) where
      fft = fft' nat
          where   fft' :: (RealFloat a) => Nat n -> RTree n (Complex a) -> RTree n (Complex a)
                  fft' Zero     = id
                  fft' (Succ _) = inDIT fft
                      where   inDIT g  = RT.toB . g . bottomSplit
  
  -- Adds the proper phase adjustments to a functor containing Complex RealFloats,
  -- and instancing Num.
  addPhase :: (Applicative f, Foldable f, LScan f, RealFloat a, Num (f (Complex a))) => f (Complex a) -> f (Complex a)
  addPhase = liftA2 (*) id phasor
    where phasor f = fst $ lproducts (pure phaseDelta)
            where phaseDelta = cis ((-pi) / fromIntegral n)
                  n          = flen f
  
  -- Gives the "length" (i.e. - number of elements in) of a Foldable.
  -- (Soon, to be provided by the Foldable class, as "length".)
  flen :: (Foldable f) => f a -> Int
  flen = foldl' (flip ((+) . const 1)) 0
  
  -- Test config.
  realData :: [[PrettyDouble]]
  realData = [  [1.0,   0.0,   0.0,   0.0]  -- Delta
              , [1.0,   1.0,   1.0,   1.0]  -- Constant
              , [1.0,  -1.0,   1.0,  -1.0]  -- Nyquist
              , [1.0,   0.0,  -1.0,   0.0]  -- Fundamental
              , [0.0,   1.0,   0.0,  -1.0]  -- Fundamental w/ 90-deg. phase lag
             ]
  complexData :: [[Complex PrettyDouble]]
  complexData = map (map (:+ 0.0)) realData
  
  myTree2 :: [a] -> RTree N2 a
  myTree2 [w, x, y, z] = RT.tree2 w x y z
  myTree2 _            = error "Something went horribly wrong!"
  
  myTree3 :: [a] -> RTree N3 a
  myTree3 [a, b, c, d, e, f, g, h] = RT.tree3 a b c d e f g h
  myTree3 _            = error "Something went horribly wrong!"
  
  myTree4 :: [a] -> RTree N4 a
  myTree4 [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = RT.tree4 a b c d e f g h i j k l m n o p
  myTree4 _            = error "Something went horribly wrong!"
  
  myTree5 :: [a] -> RTree N5 a
  myTree5 [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p,
           a', b', c', d', e', f', g', h', i', j', k', l', m', n', o', p'] =
              RT.tree5 a b c d e f g h i j k l m n o p a' b' c' d' e' f' g' h' i' j' k' l' m' n' o' p'
  myTree5 _            = error "Something went horribly wrong!"
  
  -- myTree6 :: [a] -> RTree N6 a
  -- myTree6 [a1, b1, c1, d1, e1, f1, g1, h1, i1, j1, k1, l1, m1, n1, o1, p1,
  --          a2, b2, c2, d2, e2, f2, g2, h2, i2, j2, k2, l2, m2, n2, o2, p2,
  --          a3, b3, c3, d3, e3, f3, g3, h3, i3, j3, k3, l3, m3, n3, o3, p3,
  --          a4, b4, c4, d4, e4, f4, g4, h4, i4, j4, k4, l4, m4, n4, o4, p4] =
  --             RT.tree6 a1 b1 c1 d1 e1 f1 g1 h1 i1 j1 k1 l1 m1 n1 o1 p1
  --                      a2 b2 c2 d2 e2 f2 g2 h2 i2 j2 k2 l2 m2 n2 o2 p2
  --                      a3 b3 c3 d3 e3 f3 g3 h3 i3 j3 k3 l3 m3 n3 o3 p3
  --                      a4 b4 c4 d4 e4 f4 g4 h4 i4 j4 k4 l4 m4 n4 o4 p4
  -- myTree6 _            = error "Something went horribly wrong!"
  
  -- Discrete Fourier Transform (DFT) (our "truth" reference)
  -- O(n^2)
  --
  dft :: RealFloat a => [Complex a] -> [Complex a]
  dft xs = [ sum [ x * exp((0.0 :+ (-1.0)) * 2 * pi / lenXs * fromIntegral(k * n))
                   | (x, n) <- Prelude.zip xs [0..]
                 ]
             | k <- [0..(length xs - 1)]
           ]
      where lenXs = fromIntegral $ length xs
  
  -- QuickCheck types & propositions
  newtype FFTTestVal = FFTTestVal {
      getVal :: [Complex PrettyDouble]
  } deriving (Show)
  instance Arbitrary FFTTestVal where
      arbitrary = do
          xs <- vectorOf 32 $ choose (-1.0::Double, 1.0)
          let zs = map ((:+ 0) . PrettyDouble) xs
          return $ FFTTestVal zs
  
  prop_fft_test_N2 :: FFTTestVal -> Bool
  prop_fft_test_N2 testVal = fft (myTree2 zs) == RT.fromList (dft zs)
      where zs = take 4 $ getVal testVal
  
  prop_fft_test_N3 :: FFTTestVal -> Bool
  prop_fft_test_N3 testVal = fft (myTree3 zs) == RT.fromList (dft zs)
      where zs = take 8 $ getVal testVal
  
  prop_fft_test_N4 :: FFTTestVal -> Bool
  prop_fft_test_N4 testVal = fft (myTree4 zs) == RT.fromList (dft zs)
      where zs = take 16 $ getVal testVal
  
  prop_fft_test_N5 :: FFTTestVal -> Bool
  prop_fft_test_N5 testVal = fft (myTree5 zs) == RT.fromList (dft zs)
      where zs = take 32 $ getVal testVal
  
  -- Test definitions & choice
  basicTest :: IO ()
  basicTest = forM_ complexData (\x -> do
                  putStr "\nTesting input: "
                  print x
                  putStr "Expected output: "
                  print $ dft x
                  putStr "Actual output:   "
                  print $ toList (fft (myTree2 x))
                  )
  
  -- This weirdness is required, as of GHC 7.8.
  return []
  
  runTests :: IO Bool
  runTests = $quickCheckAll
  -- End weirdness.
  
  advancedTest :: IO ()
  advancedTest = do
      allPass <- runTests -- Run QuickCheck on all prop_ functions
      unless allPass exitFailure
  
  main :: IO ()
  -- main = basicTest
  main = advancedTest
  