emd 0.1.9.0 → 0.1.10.0

raw patch · 12 files changed

+898/−254 lines, 12 filesdep +freedep +hedgehogdep +tastydep ~vector-sizedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: free, hedgehog, tasty, tasty-hedgehog, tasty-hunit

Dependency ranges changed: vector-sized

API changes (from Hackage documentation)

Files

• CHANGELOG.md +19/−0
• emd.cabal +22/−6
• src/Numeric/EMD.hs +18/−219
• src/Numeric/EMD/Internal/Pipe.hs +224/−0
• src/Numeric/EMD/Internal/Sift.hs +384/−0
• src/Numeric/EMD/Internal/Spline.hs +1/−1
• src/Numeric/HHT.hs +2/−2
• test/Spec.hs +8/−26
• test/Tests/EMD.hs +82/−0
• test/Tests/HHT.hs +31/−0
• test/Tests/Spline.hs +28/−0
• test/Tests/Util.hs +79/−0

@@ -1,6 +1,25 @@ Changelog ========= +Version 0.1.10.0+---------------++*October 19, 2019*++<https://github.com/mstksg/emd/releases/tag/v0.1.10.0>++*   Sifting condition checking system has been revamped to be more flexible.+*   New sift conditions added:+    *   Projections based on functions:+        *   Energy Difference function (for Energy Difference Tracking)+        *   Mean envelope RMS+    *   For the above, we can now stop sifting based on:+        *   Those projections reaching a certain value+        *   Successive values of those projections reaching a certain squared+            difference (Cauchy-inspired convergence)+    *   S-Number Criterion+*   Add test suites+ Version 0.1.9.0 --------------- 

@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 84c1c690a80a0ace06a7280dfac722ed522cdd34fb042042999e67d7a1217442+-- hash: 5230a7607e7d3f26121daa88b982f8f58087a236adc506e45eadf1a038581d7a  name:           emd-version:        0.1.9.0+version:        0.1.10.0 synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform description:    Empirical Mode decomposition and Hilbert-Huang Transform in pure                 Haskell.@@ -36,12 +36,14 @@       Numeric.EMD.Internal.Spline       Numeric.HHT   other-modules:-      Numeric.EMD.Internal.Tridiagonal       Numeric.EMD.Internal.Extrema+      Numeric.EMD.Internal.Pipe+      Numeric.EMD.Internal.Sift+      Numeric.EMD.Internal.Tridiagonal       Numeric.HHT.Internal.FFT   hs-source-dirs:       src-  ghc-options: -Wall -Wredundant-constraints -Wcompat+  ghc-options: -Wall -Wredundant-constraints -Wcompat -Werror=incomplete-patterns   build-depends:       array     , base >=4.11 && <5@@ -52,6 +54,7 @@     , deepseq     , fft     , finite-typelits+    , free     , ghc-typelits-knownnat     , ghc-typelits-natnormalise     , transformers@@ -64,15 +67,28 @@   type: exitcode-stdio-1.0   main-is: Spec.hs   other-modules:+      Tests.EMD+      Tests.HHT+      Tests.Spline+      Tests.Util       Paths_emd   hs-source-dirs:       test-  ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N+  ghc-options: -Wall -Wredundant-constraints -Wcompat -Werror=incomplete-patterns -threaded -rtsopts -with-rtsopts=-N   build-depends:       HUnit     , base >=4.11 && <5     , containers     , emd+    , ghc-typelits-knownnat+    , ghc-typelits-natnormalise+    , hedgehog+    , statistics+    , tasty+    , tasty-hedgehog+    , tasty-hunit+    , typelits-witnesses+    , vector-sized   default-language: Haskell2010  benchmark emd-bench@@ -82,7 +98,7 @@       Paths_emd   hs-source-dirs:       bench-  ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N -O2+  ghc-options: -Wall -Wredundant-constraints -Wcompat -Werror=incomplete-patterns -threaded -rtsopts -with-rtsopts=-N -O2   build-depends:       base >=4.11 && <5     , criterion

@@ -2,6 +2,7 @@ {-# LANGUAGE DeriveGeneric                            #-} {-# LANGUAGE GADTs                                    #-} {-# LANGUAGE LambdaCase                               #-}+{-# LANGUAGE RankNTypes                               #-} {-# LANGUAGE RecordWildCards                          #-} {-# LANGUAGE ScopedTypeVariables                      #-} {-# LANGUAGE TypeApplications                         #-}@@ -33,9 +34,9 @@ -- 2.  We provide a vector of size of at least one. -- -- There are many functions to convert unsized vectors to sized vectors in--- "Data.Vector.Sized" and associated modules, including 'toSized' (for--- when you know the size at compile-time) and 'withSized' (for when you--- don't).+-- "Data.Vector.Sized" and associated modules, including+-- 'Data.Vector.Sized.toSized' (for when you know the size at compile-time)+-- and 'Data.Vector.Sized.withSized' (for when you don't). --  module Numeric.EMD (@@ -45,116 +46,28 @@   , emd'   , iemd   , EMD(..)-  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)+  -- ** Configuration+  , EMDOpts(..), defaultEO+  , BoundaryHandler(..)+  , SiftCondition(..), SiftProjection(..), defaultSC+  , scEnergyDiff+  , SplineEnd(..)   -- * Internal   , sift, SiftResult(..)   , envelopes   ) where  import           Control.DeepSeq-import           Control.Monad import           Control.Monad.IO.Class-import           Data.Default.Class-import           Data.Finite import           Data.Functor.Identity import           Data.List-import           GHC.Generics                 (Generic)+import           GHC.Generics              (Generic) import           GHC.TypeNats-import           Numeric.EMD.Internal.Extrema-import           Numeric.EMD.Internal.Spline+import           Numeric.EMD.Internal.Sift import           Text.Printf-import qualified Data.Binary                  as Bi-import qualified Data.Map                     as M-import qualified Data.Vector.Generic          as VG-import qualified Data.Vector.Generic.Sized    as SVG---- | Options for EMD composition.-data EMDOpts a = EO { eoSiftCondition   :: SiftCondition a  -- ^ stop condition for sifting-                    , eoSplineEnd       :: SplineEnd a      -- ^ end conditions for envelope splines-                    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary-                    }-  deriving (Show, Eq, Ord, Generic)--data BoundaryHandler-    -- | Clamp envelope at end points (Matlab implementation)-    = BHClamp-    -- | Extend boundaries symmetrically-    | BHSymmetric-  deriving (Show, Eq, Ord, Generic)--    -- -- | Extend boundaries assuming global periodicity-    -- -- | BHPeriodic---- | @since 0.1.3.0-instance Bi.Binary BoundaryHandler---- | @since 0.1.3.0-instance Bi.Binary a => Bi.Binary (EMDOpts a)---- | Default 'EMDOpts'-defaultEO :: Fractional a => EMDOpts a-defaultEO = EO { eoSiftCondition   = defaultSC-               , eoSplineEnd       = SENatural-               , eoBoundaryHandler = Just BHSymmetric-               }---- | @since 0.1.3.0-instance Fractional a => Default (EMDOpts a) where-    def = defaultEO---- | Stop conditions for sifting process------ Data type is lazy in its fields, so this infinite data type:------ @--- nTimes n = SCTimes n `SCOr` nTimes (n + 1)--- @------ will be treated identically as:------ @--- nTimes = SCTimes--- @-data SiftCondition a-    -- | Stop using standard SD method-    = SCStdDev !a-    -- | Stop after a fixed number of sifting iterations-    | SCTimes !Int-    -- | One or the other-    | SCOr (SiftCondition a) (SiftCondition a)-    -- | Stop when both conditions are met-    | SCAnd (SiftCondition a) (SiftCondition a)-  deriving (Show, Eq, Ord, Generic)---- | @since 0.1.3.0-instance Bi.Binary a => Bi.Binary (SiftCondition a)---- | @since 0.1.3.0-instance Fractional a => Default (SiftCondition a) where-    def = defaultSC---- | Default 'SiftCondition'-defaultSC :: Fractional a => SiftCondition a-defaultSC = SCStdDev 0.3 `SCOr` SCTimes 50     -- R package uses SCTimes 20, Matlab uses no limit--- defaultSC = SCStdDev 0.3---- | 'True' if stop-testCondition-    :: (VG.Vector v a, Fractional a, Ord a)-    => SiftCondition a-    -> Int-    -> SVG.Vector v n a-    -> SVG.Vector v n a-    -> Bool-testCondition tc i v v' = go tc-  where-    sd = SVG.sum $ SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v'-    go = \case-      SCStdDev t -> sd <= t-      SCTimes l  -> i >= l-      SCOr  f g  -> go f || go g-      SCAnd f g  -> go f && go g-    eps = 0.0000001+import qualified Data.Binary               as Bi+import qualified Data.Vector.Generic       as VG+import qualified Data.Vector.Generic.Sized as SVG  -- | An @'EMD' v n a@ is an Empirical Mode Decomposition of a time series -- with @n@ items of type @a@ stored in a vector @v@.@@ -186,7 +99,7 @@ -- 1.  The resulting 'EMD' contains IMFs that are all the same length as --     the input vector -- 2.  We provide a vector of size of at least one.-emd :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)+emd :: (VG.Vector v a, KnownNat n, Floating a, Ord a)     => EMDOpts a     -> SVG.Vector v (n + 1) a     -> EMD v (n + 1) a@@ -195,7 +108,7 @@ -- | 'emd', but tracing results to stdout as IMFs are found.  Useful for -- debugging to see how long each step is taking. emdTrace-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a, MonadIO m)+    :: (VG.Vector v a, KnownNat n, Floating a, Ord a, MonadIO m)     => EMDOpts a     -> SVG.Vector v (n + 1) a     -> m (EMD v (n + 1) a)@@ -205,7 +118,7 @@  -- | 'emd' with a callback for each found IMF. emd'-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a, Applicative m)+    :: (VG.Vector v a, KnownNat n, Floating a, Ord a, Applicative m)     => (SiftResult v (n + 1) a -> m r)     -> EMDOpts a     -> SVG.Vector v (n + 1) a@@ -228,117 +141,3 @@     => EMD v n a     -> SVG.Vector v n a iemd EMD{..} = foldl' (SVG.zipWith (+)) emdResidual emdIMFs---- | The result of a sifting operation.  Each sift either yields--- a residual, or a new IMF.-data SiftResult v n a = SRResidual !(SVG.Vector v n a)-                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of sifting iterations---- | Iterated sifting process, used to produce either an IMF or a residual.-sift-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => EMDOpts a-    -> SVG.Vector v (n + 1) a-    -> SiftResult v (n + 1) a-sift EO{..} = go 1-  where-    go !i !v = case sift' eoSplineEnd eoBoundaryHandler v of-      Nothing -> SRResidual v-      Just !v'-        | testCondition eoSiftCondition i v v' -> SRIMF v' i-        | otherwise                            -> go (i + 1) v'---- | Single sift-sift'-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd a-    -> Maybe BoundaryHandler-    -> SVG.Vector v (n + 1) a-    -> Maybe (SVG.Vector v (n + 1) a)-sift' se bh v = go <$> envelopes se bh v-  where-    go (mins, maxs) = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs---- | Returns cubic splines of local minimums and maximums.  Returns--- 'Nothing' if there are not enough local minimum or maximums to create--- the splines.-envelopes-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd a-    -> Maybe BoundaryHandler-    -> SVG.Vector v (n + 1) a-    -> Maybe (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a)-envelopes se bh xs = do-    when (bh == Just BHClamp) $ do-      guard (M.size mins > 1)-      guard (M.size maxs > 1)-    (,) <$> splineAgainst se emin mins-        <*> splineAgainst se emax maxs-  where-    -- minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]-    (mins,maxs) = extrema xs-    (emin,emax) = case bh of-      Nothing  -> mempty-      Just bh' -> extendExtrema xs bh' (mins,maxs)-    --   | isJust bh = (mins `M.union` minMax, maxs `M.union` minMax)-    --   | otherwise = (mins, maxs)--extendExtrema-    :: forall v n a. (VG.Vector v a, KnownNat n)-    => SVG.Vector v (n + 1) a-    -> BoundaryHandler-    -> (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)-    -> (M.Map Int a, M.Map Int a)-    -- (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)-extendExtrema xs = \case-    BHClamp     -> const (firstLast, firstLast)-    BHSymmetric -> \(mins, maxs) ->-      let addFirst = case (flippedMin, flippedMax) of-              (Nothing      , Nothing      ) -> mempty-              -- first point is local maximum-              (Just (_,mn)  , Nothing      ) -> (mn        , firstPoint)-              -- first point is local minimum-              (Nothing      , Just (_,mx)  ) -> (firstPoint, mx        )-              (Just (mni,mn), Just (mxi,mx))-                | mni < mxi                  -> (mn        , firstPoint)-                | otherwise                  -> (firstPoint, mx        )-            where-              flippedMin = flip fmap (M.lookupMin mins) $ \(minIx, minVal) ->-                (minIx, M.singleton (negate (fromIntegral minIx)) minVal)-              flippedMax = flip fmap (M.lookupMin maxs) $ \(maxIx, maxVal) ->-                (maxIx, M.singleton (negate (fromIntegral maxIx)) maxVal)-          addLast = case (flippedMin, flippedMax) of-              (Nothing      , Nothing      ) -> mempty-              -- last point is local maximum-              (Just (_,mn)  , Nothing      ) -> (mn        , lastPoint )-              -- last point is local minimum-              (Nothing      , Just (_,mx)  ) -> (lastPoint , mx        )-              (Just (mni,mn), Just (mxi,mx))-                | mni > mxi                  -> (mn        , lastPoint )-                | otherwise                  -> (lastPoint , mx        )-            where-              flippedMin = flip fmap (M.lookupMax mins) $ \(minIx, minVal) ->-                (minIx, M.singleton (extendSym (fromIntegral minIx)) minVal)-              flippedMax = flip fmap (M.lookupMax maxs) $ \(maxIx, maxVal) ->-                (maxIx, M.singleton (extendSym (fromIntegral maxIx)) maxVal)-      in  addFirst `mappend` addLast-  where-    lastIx = fromIntegral $ maxBound @(Finite n)-    firstPoint = M.singleton 0 (SVG.head xs)-    lastPoint  = M.singleton lastIx (SVG.last xs)-    firstLast  = firstPoint `mappend` lastPoint-    extendSym i = 2 * lastIx - i---- | Build a splined vector against a map of control points.-splineAgainst-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd a-    -> M.Map Int a              -- ^ extensions-    -> M.Map (Finite n) a-    -> Maybe (SVG.Vector v n a)-splineAgainst se ext = fmap go-                     . makeSpline se-                     . mappend (M.mapKeysMonotonic fromIntegral ext)-                     . M.mapKeysMonotonic fromIntegral-  where-    go spline = SVG.generate (sampleSpline spline . fromIntegral)

@@ -0,0 +1,224 @@+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TemplateHaskell            #-}++module Numeric.EMD.Internal.Pipe (+    Pipe+  , (.|)+  , runPipe+  , awaitEither, await, awaitSurely+  , repeatM, unfoldP, unfoldPForever, iterateP, sourceList+  , awaitForever, mapP, mapMP+  , dropP+  , foldrP, sinkList+  , ZipSink(..)+  ) where++import           Control.Applicative+import           Control.Monad+import           Control.Monad.Free.Class+import           Control.Monad.Free.TH+import           Control.Monad.Trans.Class+import           Control.Monad.Trans.Free        (FreeT(..), FreeF(..))+import           Control.Monad.Trans.Free.Church+import           Data.Foldable+import           Data.Void++data PipeF i o u a =+      PAwaitF (i -> a) (u -> a)+    | PYieldF o a+  deriving Functor++makeFree ''PipeF++-- | Similar to Conduit+--+-- *  @i@: Type of input stream+-- *  @o@: Type of output stream+-- *  @u@: Type of the /result/ of the upstream pipe (Outputted when+--    upstream pipe finishes)+-- *  @m@: Underlying monad+-- *  @a@: Result type (Outputted when finished)+--+-- Some specializations:+--+-- *  A pipe is a /producer/ if @i@ is @()@: it doesn't need anything to go+--    pump out items.+--+--    If a pipe is producer and @a@ is 'Void', it means that it will+--    produce infinitely.+--+-- *  A pipe is a /consumer/ if @o@ is 'Void': it will never yield anything+--    else downstream.+--+-- *  Normally you can ask for input upstream with 'await', which returns+--    'Nothing' if the pipe upstream stops producing.  However, if @u@ is+--    'Void', it means that the pipe upstream will never stop, so you can+--    use 'awaitSurely' to get a guaranteed answer.+newtype Pipe i o u m a = Pipe { pipeFree :: FT (PipeF i o u) m a }+  deriving (Functor, Applicative, Monad, MonadTrans, MonadFree (PipeF i o u))++awaitEither :: Pipe i o u m (Either i u)+awaitEither = pAwaitF++yield :: o -> Pipe i o u m ()+yield = pYieldF++await :: Pipe i o u m (Maybe i)+await = either Just (const Nothing) <$> awaitEither++awaitSurely :: Pipe i o Void m i+awaitSurely = either id absurd <$> awaitEither++runPipe :: forall u m a. Monad m => Pipe () Void u m a -> m a+runPipe = iterT go . pipeFree+  where+    go :: PipeF () Void u (m a) -> m a+    go = \case+      PAwaitF f _ -> f ()+      PYieldF o _ -> absurd o++-- can this be done without going through FreeT?+(.|) :: forall a b c u m v r. Monad m => Pipe a b u m v -> Pipe b c v m r -> Pipe a c u m r+Pipe p .| Pipe q = Pipe $ toFT $ compPipe_ (fromFT p) (fromFT q)++compPipe_+    :: forall a b c u v m r. (Monad m)+    => FreeT (PipeF a b u) m v+    -> FreeT (PipeF b c v) m r+    -> FreeT (PipeF a c u) m r+compPipe_ p q = FreeT $ runFreeT q >>= \case+    Pure x             -> pure . Pure $ x+    Free (PAwaitF f g) -> runFreeT p >>= \case+      Pure x'              -> runFreeT $ compPipe_ p  (g x')+      Free (PAwaitF f' g') -> pure . Free $ PAwaitF ((`compPipe_` q) . f')+                                                    ((`compPipe_` q) . g')+      Free (PYieldF x' y') -> runFreeT $ compPipe_ y' (f x')+    Free (PYieldF x y) -> pure . Free $ PYieldF x (compPipe_ p y)+infixr 2 .|++unfoldP :: (b -> Maybe (a, b)) -> b -> Pipe i a u m ()+unfoldP f = go+  where+    go z = case f z of+      Nothing      -> pure ()+      Just (x, z') -> yield x *> go z'++unfoldPForever :: (b -> (a, b)) -> b -> Pipe i a u m r+unfoldPForever f = go+  where+    go z = yield x *> go z'+      where+        (x, z') = f z++iterateP :: (a -> a) -> a -> Pipe i a u m r+iterateP f = unfoldPForever (join (,) . f)++sourceList :: Foldable t => t a -> Pipe i a u m ()+sourceList = traverse_ yield++repeatM :: Monad m => m o -> Pipe i o u m u+repeatM x = go+  where+    go = (yield =<< lift x) *> go++awaitForever :: (i -> Pipe i o u m a) -> Pipe i o u m u+awaitForever f = go+  where+    go = awaitEither >>= \case+      Left x  -> f x *> go+      Right x -> pure x++-- finishPipe+--     :: u+--     -> Pipe i o u    m a+--     -> Pipe i o Void m a++mapP :: (a -> b) -> Pipe a b u m u+mapP f = awaitForever (yield . f)++mapMP :: Monad m => (a -> m b) -> Pipe a b u m u+mapMP f = awaitForever ((yield =<<) . lift . f)++dropP :: Int -> Pipe i o u m ()+dropP n = replicateM_ n await++foldrP :: (a -> b -> b) -> b -> Pipe a Void u m b+foldrP f z = go+  where+    go = await >>= \case+      Nothing -> pure z+      Just x  -> f x <$> go++sinkList :: Pipe i Void u m [i]+sinkList = foldrP (:) []++newtype ZipSink i u m a = ZipSink { getZipSink :: Pipe i Void u m a }+  deriving Functor++zipSink_+    :: Monad m+    => FreeT (PipeF i Void u) m (a -> b)+    -> FreeT (PipeF i Void u) m a+    -> FreeT (PipeF i Void u) m b+zipSink_ p q = FreeT $ go <$> runFreeT p <*> runFreeT q+  where+    go = \case+      Pure x             -> \case+        Pure x'              -> Pure $ x x'+        Free (PAwaitF f' g') -> Free $ PAwaitF (zipSink_ p . f') (zipSink_ p . g')+        Free (PYieldF x' _ ) -> absurd x'+      Free (PAwaitF f g) -> \case+        Pure _               -> Free $ PAwaitF ((`zipSink_` q) . f) ((`zipSink_` q) . g)+        Free (PAwaitF f' g') -> Free $ PAwaitF (zipSink_ <$> f <*> f') (zipSink_ <$> g <*> g')+        Free (PYieldF x' _ ) -> absurd x'+      Free (PYieldF x _) -> absurd x++altSink_+    :: Monad m+    => FreeT (PipeF i Void u) m a+    -> FreeT (PipeF i Void u) m a+    -> FreeT (PipeF i Void u) m a+altSink_ p q = FreeT $ go <$> runFreeT p <*> runFreeT q+  where+    go = \case+      Pure x             -> \_ -> Pure x+      Free (PAwaitF f g) -> \case+        Pure x'              -> Pure x'+        Free (PAwaitF f' g') -> Free $ PAwaitF (altSink_ <$> f <*> f') (altSink_ <$> g <*> g')+        Free (PYieldF x' _ ) -> absurd x'+      Free (PYieldF x _) -> absurd x++zipSink+    :: Monad m+    => Pipe i Void u m (a -> b)+    -> Pipe i Void u m a+    -> Pipe i Void u m b+zipSink (Pipe p) (Pipe q) = Pipe $ toFT $ zipSink_ (fromFT p) (fromFT q)++altSink+    :: Monad m+    => Pipe i Void u m a+    -> Pipe i Void u m a+    -> Pipe i Void u m a+altSink (Pipe p) (Pipe q) = Pipe $ toFT $ altSink_ (fromFT p) (fromFT q)++-- | '<*>' = distribute input to all, and return result when they finish+--+-- 'pure' = immediately finish+instance Monad m => Applicative (ZipSink i u m) where+    pure = ZipSink . pure+    ZipSink p <*> ZipSink q = ZipSink $ zipSink p q++-- | '<|>' = distribute input to all, and return the first result that+-- finishes+--+-- 'empty' = never finish+instance Monad m => Alternative (ZipSink i u m) where+    empty = ZipSink go+      where+        go = forever await+    ZipSink p <|> ZipSink q = ZipSink $ altSink p q

@@ -0,0 +1,384 @@+{-# LANGUAGE BangPatterns                             #-}+{-# LANGUAGE DeriveGeneric                            #-}+{-# LANGUAGE GADTs                                    #-}+{-# LANGUAGE LambdaCase                               #-}+{-# LANGUAGE RankNTypes                               #-}+{-# LANGUAGE RecordWildCards                          #-}+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeInType                               #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}+++module Numeric.EMD.Internal.Sift (+    EMDOpts(..), defaultEO+  , BoundaryHandler(..)+  , SiftCondition(..), SiftProjection(..), defaultSC+  , scEnergyDiff+  , SplineEnd(..)+  -- * Internal+  , sift, SiftResult(..)+  , envelopes+  ) where++import           Control.Applicative+import           Control.Monad+import           Control.Monad.Trans.State+import           Data.Default.Class+import           Data.Finite+import           Data.Void+import           GHC.Generics                 (Generic)+import           GHC.TypeNats+import           Numeric.EMD.Internal.Extrema+import           Numeric.EMD.Internal.Pipe+import           Numeric.EMD.Internal.Spline+import qualified Data.Binary                  as Bi+import qualified Data.Map                     as M+import qualified Data.Vector.Generic          as VG+import qualified Data.Vector.Generic.Sized    as SVG++-- | Options for EMD composition.+data EMDOpts a = EO { eoSiftCondition   :: SiftCondition a  -- ^ stop condition for sifting+                    , eoSplineEnd       :: SplineEnd a      -- ^ end conditions for envelope splines+                    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary+                    }+  deriving (Show, Eq, Ord, Generic)++-- | @since 0.1.3.0+instance Bi.Binary a => Bi.Binary (EMDOpts a)++-- | Default 'EMDOpts'+defaultEO :: Fractional a => EMDOpts a+defaultEO = EO { eoSiftCondition   = defaultSC+               , eoSplineEnd       = SENatural+               , eoBoundaryHandler = Just BHSymmetric+               }++-- | @since 0.1.3.0+instance Fractional a => Default (EMDOpts a) where+    def = defaultEO+++-- | Boundary conditions for splines.+data BoundaryHandler+    -- | Clamp envelope at end points (Matlab implementation)+    = BHClamp+    -- | Extend boundaries symmetrically+    | BHSymmetric+  deriving (Show, Eq, Ord, Generic)++-- | @since 0.1.3.0+instance Bi.Binary BoundaryHandler++-- | Stop conditions for sifting process+--+-- Data type is lazy in its fields, so this infinite data type:+--+-- @+-- nTimes n = SCTimes n `SCOr` nTimes (n + 1)+-- @+--+-- will be treated identically as:+--+-- @+-- nTimes = SCTimes+-- @+data SiftCondition a+    -- | Stop using standard SD method+    = SCStdDev !a+    -- | When the difference between successive items reaches a given threshold+    -- \(\tau\)+    --+    -- \[+    -- \frac{\left(f(t-1) - f(t)\right)^2}{f^2(t-1)} < \tau+    -- \]+    --+    -- @since 0.1.10.0+    | SCCauchy SiftProjection !a+    -- | When the value reaches a given threshold \(\tau\)+    --+    -- \[+    -- f(t) < \tau+    -- \]+    --+    -- @since 0.1.10.0+    | SCProj   SiftProjection !a+    -- | S-condition criteria.+    --+    -- The S-number is the length of current streak where number of extrema+    -- or zero crossings all differ at most by one.+    --+    -- Stop sifting when the S-number reaches a given amount.+    --+    -- @since 0.1.10.0+    | SCSCond !Int+    -- | Stop after a fixed number of sifting iterations+    | SCTimes !Int+    -- | One or the other+    | SCOr (SiftCondition a) (SiftCondition a)+    -- | Stop when both conditions are met+    | SCAnd (SiftCondition a) (SiftCondition a)+  deriving (Show, Eq, Ord, Generic)++-- | A projection of sifting data.  Used as a part of 'SiftCondition' to+-- describe 'SCCauchy' and 'SCProj'.+--+-- @since 0.1.10.0+data SiftProjection+    -- | The root mean square of the envelope means+    = SPEnvMeanSum+    -- | The "energy difference" quotient (Cheng, Yu, Yang 2005)+    | SPEnergyDiff+  deriving (Show, Eq, Ord, Generic)++instance Bi.Binary SiftProjection++-- | @since 0.1.3.0+instance Bi.Binary a => Bi.Binary (SiftCondition a)++-- | @since 0.1.3.0+instance Fractional a => Default (SiftCondition a) where+    def = defaultSC++-- | Default 'SiftCondition'+defaultSC :: Fractional a => SiftCondition a+defaultSC = SCStdDev 0.3 `SCOr` SCTimes 50     -- R package uses SCTimes 20, Matlab uses no limit+++-- | Cheng, Yu, Yang suggest pairing together an energy difference+-- threshold with a threshold for mean envelope RMS.  This is a convenience+-- function to construct that pairing.+scEnergyDiff+    :: a                -- ^ Threshold for Energy Difference+    -> a                -- ^ Threshold for mean envelope RMS+    -> SiftCondition a+scEnergyDiff s t = SCProj SPEnergyDiff s `SCAnd` SCProj SPEnvMeanSum t+++-- | The result of a sifting operation.  Each sift either yields+-- a residual, or a new IMF.+data SiftResult v n a = SRResidual !(SVG.Vector v n a)+                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of sifting iterations++-- | Result of a single sift+data SingleSift v n a = SingleSift+    { ssRes    :: !(SVG.Vector v n a)+    , ssMinEnv :: !(SVG.Vector v n a)+    , ssMaxEnv :: !(SVG.Vector v n a)+    }++type Sifter v n m a = Pipe (SingleSift v n a) Void Void m ()++siftTimes :: Int -> Sifter v n m a+siftTimes n = dropP (n - 1) >> void awaitSurely++siftProj :: (SingleSift v n a -> Bool) -> Sifter v n m a+siftProj p = go+  where+    go = do+      v <- awaitSurely+      unless (p v) go++siftPairs :: (SingleSift v n a -> SingleSift v n a -> Bool) -> Sifter v n m a+siftPairs p = go =<< awaitSurely+  where+    go s = do+      s' <- awaitSurely+      unless (p s s') (go s')++siftStdDev :: forall v n m a. (VG.Vector v a, Fractional a, Ord a) => a -> Sifter v n m a+siftStdDev t = siftPairs $ \(SingleSift v _ _) (SingleSift v' _ _) ->+    SVG.sum (SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v')+      <= t+  where+    eps = 0.0000001++siftCauchy+    :: (Fractional b, Ord b)+    => (SingleSift v n a -> b)+    -> b+    -> Sifter v n m a+siftCauchy p t = siftPairs $ \s s' ->+  let ps  = p s+      ps' = p s'+      δ   = ps' - ps+  in  ((δ * δ) / (ps * ps)) <= t++siftSCond :: (VG.Vector v a, KnownNat n, Fractional a, Ord a) => Int -> Sifter v (n + 1) m a+siftSCond n = go []+  where+    go cxs = do+      v <- awaitSurely+      let cx   = crossCount $ ssRes v+          done = all ((<= 1) . abs . subtract cx) cxs+      unless done $+        go (take (n - 1) (cx : cxs))+    crossCount xs = M.size mins + M.size maxs + crosses+      where+        (mins, maxs) = extrema xs+        crosses = fst . flip execState (0, Nothing) . flip SVG.mapM_ xs $ \x -> modify $ \(!i, !y) ->+          let xPos = x > 0+              i'   = case y of+                       Nothing -> i+                       Just y'+                         | xPos == y' -> i+                         | otherwise  -> i + 1+          in  (i', Just xPos)++siftOr :: Monad m => Sifter v n m a -> Sifter v n m a -> Sifter v n m a+siftOr p q = getZipSink $ ZipSink p <|> ZipSink q++siftAnd :: Monad m => Sifter v n m a -> Sifter v n m a -> Sifter v n m a+siftAnd p q = getZipSink $ ZipSink p *> ZipSink q++toSifter+    :: (VG.Vector v a, KnownNat n, Monad m, Floating a, Ord a)+    => SVG.Vector v (n + 1) a+    -> SiftCondition a+    -> Sifter v (n + 1) m a+toSifter v0 = go+  where+    go = \case+      SCStdDev x -> siftStdDev x+      SCCauchy p x -> siftCauchy (toProj p v0) x+      SCProj   p x -> siftProj ((<= x) . toProj p v0)+      SCSCond  n   -> siftSCond n+      SCTimes  i -> siftTimes i+      SCOr p q   -> siftOr (go p) (go q)+      SCAnd p q  -> siftAnd (go p) (go q)++toProj+    :: (VG.Vector v a, Floating a)+    => SiftProjection+    -> SVG.Vector v n a+    -> SingleSift v n a+    -> a+toProj = \case+    SPEnvMeanSum -> \_ SingleSift{..} ->+      sqrt . squareMag $ SVG.zipWith (\x y -> (x + y) / 2) ssMinEnv ssMaxEnv+    SPEnergyDiff -> \v0 ->+      let eX = squareMag v0+      in  \SingleSift{..} ->+            let eTot = squareMag ssRes - squareMag (SVG.zipWith (-) v0 ssRes)+            in  abs $ eX - eTot+  where+    squareMag = SVG.foldl' (\s x -> s + x*x) 0+++-- | Iterated sifting process, used to produce either an IMF or a residual.+sift+    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a, Ord a)+    => EMDOpts a+    -> SVG.Vector v (n + 1) a+    -> SiftResult v (n + 1) a+sift EO{..} v0 = case execStateT (runPipe sifterPipe) (0, v0) of+    Left  v        -> SRResidual v+    Right (!i, !v) -> SRIMF v i+  where+    sifterPipe = repeatM go+              .| toSifter v0 eoSiftCondition+    go = StateT $ \(!i, !v) ->+      case sift' eoSplineEnd eoBoundaryHandler v of+        Nothing                -> Left v+        Just ss@SingleSift{..} -> Right (ss, (i + 1, ssRes))++-- | Single sift+sift'+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)+    => SplineEnd a+    -> Maybe BoundaryHandler+    -> SVG.Vector v (n + 1) a+    -> Maybe (SingleSift v (n + 1) a)+sift' se bh v = do+    (mins, maxs) <- envelopes se bh v+    pure SingleSift+      { ssRes    = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs+      , ssMinEnv = mins+      , ssMaxEnv = maxs+      }++-- | Returns cubic splines of local minimums and maximums.  Returns+-- 'Nothing' if there are not enough local minimum or maximums to create+-- the splines.+envelopes+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)+    => SplineEnd a+    -> Maybe BoundaryHandler+    -> SVG.Vector v (n + 1) a+    -> Maybe (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a)+envelopes se bh xs = do+    when (bh == Just BHClamp) $ do+      guard (M.size mins > 1)+      guard (M.size maxs > 1)+    (,) <$> splineAgainst se emin mins+        <*> splineAgainst se emax maxs+  where+    -- minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]+    (mins,maxs) = extrema xs+    (emin,emax) = case bh of+      Nothing  -> mempty+      Just bh' -> extendExtrema xs bh' (mins,maxs)+    --   | isJust bh = (mins `M.union` minMax, maxs `M.union` minMax)+    --   | otherwise = (mins, maxs)++extendExtrema+    :: forall v n a. (VG.Vector v a, KnownNat n)+    => SVG.Vector v (n + 1) a+    -> BoundaryHandler+    -> (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)+    -> (M.Map Int a, M.Map Int a)+    -- (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)+extendExtrema xs = \case+    BHClamp     -> const (firstLast, firstLast)+    BHSymmetric -> \(mins, maxs) ->+      let addFirst = case (flippedMin, flippedMax) of+              (Nothing      , Nothing      ) -> mempty+              -- first point is local maximum+              (Just (_,mn)  , Nothing      ) -> (mn        , firstPoint)+              -- first point is local minimum+              (Nothing      , Just (_,mx)  ) -> (firstPoint, mx        )+              (Just (mni,mn), Just (mxi,mx))+                | mni < mxi                  -> (mn        , firstPoint)+                | otherwise                  -> (firstPoint, mx        )+            where+              flippedMin = flip fmap (M.lookupMin mins) $ \(minIx, minVal) ->+                (minIx, M.singleton (negate (fromIntegral minIx)) minVal)+              flippedMax = flip fmap (M.lookupMin maxs) $ \(maxIx, maxVal) ->+                (maxIx, M.singleton (negate (fromIntegral maxIx)) maxVal)+          addLast = case (flippedMin, flippedMax) of+              (Nothing      , Nothing      ) -> mempty+              -- last point is local maximum+              (Just (_,mn)  , Nothing      ) -> (mn        , lastPoint )+              -- last point is local minimum+              (Nothing      , Just (_,mx)  ) -> (lastPoint , mx        )+              (Just (mni,mn), Just (mxi,mx))+                | mni > mxi                  -> (mn        , lastPoint )+                | otherwise                  -> (lastPoint , mx        )+            where+              flippedMin = flip fmap (M.lookupMax mins) $ \(minIx, minVal) ->+                (minIx, M.singleton (extendSym (fromIntegral minIx)) minVal)+              flippedMax = flip fmap (M.lookupMax maxs) $ \(maxIx, maxVal) ->+                (maxIx, M.singleton (extendSym (fromIntegral maxIx)) maxVal)+      in  addFirst `mappend` addLast+  where+    lastIx = fromIntegral $ maxBound @(Finite n)+    firstPoint = M.singleton 0 (SVG.head xs)+    lastPoint  = M.singleton lastIx (SVG.last xs)+    firstLast  = firstPoint `mappend` lastPoint+    extendSym i = 2 * lastIx - i++-- | Build a splined vector against a map of control points.+splineAgainst+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)+    => SplineEnd a+    -> M.Map Int a              -- ^ extensions+    -> M.Map (Finite n) a+    -> Maybe (SVG.Vector v n a)+splineAgainst se ext = fmap go+                     . makeSpline se+                     . mappend (M.mapKeysMonotonic fromIntegral ext)+                     . M.mapKeysMonotonic fromIntegral+  where+    go spline = SVG.generate (sampleSpline spline . fromIntegral)

@@ -1,6 +1,6 @@ {-# LANGUAGE ApplicativeDo                            #-} {-# LANGUAGE DataKinds                                #-}-{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DeriveGeneric                            #-} {-# LANGUAGE GADTs                                    #-} {-# LANGUAGE RecordWildCards                          #-} {-# LANGUAGE ScopedTypeVariables                      #-}

@@ -204,8 +204,8 @@ -- Takes a "binning" function to allow you to specify how specific you want -- your frequencies to be. ----- See 'hhtSparseSpetrum' for a sparser version, and 'hhtDenseSpectrum' for--- a denser version.+-- See 'hhtSparseSpectrum' for a sparser version, and 'hhtDenseSpectrum'+-- for a denser version. hhtSpectrum     :: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Num a)     => (a -> k)     -- ^ binning function.  takes rev/tick freq between 0 and 1.

@@ -1,30 +1,12 @@-{-# LANGUAGE TypeApplications #-} -import           Control.Monad-import           Data.Maybe-import           Numeric.EMD.Internal.Spline-import           Test.HUnit-import qualified Data.Map                as M-import qualified Data.Set                as S+import           Test.Tasty+import           Tests.EMD+import           Tests.HHT+import           Tests.Spline  main :: IO ()-main = void . runTestTT $ TestList-    [ "Sine spline" ~: splineTest+main = defaultMain $ testGroup "Tests"+    [ splineTest+    , emdTests+    , hhtTests     ]--splineTest :: Assertion-splineTest = do-    expected <- map (read @Double) . lines <$> readFile "test-data/sintest.csv"-    roundOut expected @=? roundOut samples-  where-    roundOut :: [Double] -> [Double]-    roundOut = map $ (/ 10e12) . fromInteger . round . (* 10e12)-    spline :: Spline Double-    spline = fromJust-           . makeSpline SENotAKnot-           . M.fromSet sin-           . S.fromList-           $ [0, 1, 2.5, 3.6, 5, 7, 8.1, 10]-    samples :: [Double]-    samples = sampleSpline spline . (/ 4) . fromInteger <$> [0..40]-

@@ -0,0 +1,82 @@+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TemplateHaskell                          #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeInType                               #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}++module Tests.EMD (+    emdTests+  ) where++import           Control.Monad+import           Data.Functor.Identity+import           Data.Proxy+import           GHC.TypeNats+import           Hedgehog+import           Numeric.EMD+import           Test.Tasty+import           Tests.Util+import qualified Hedgehog.Range        as Range++emdTests :: TestTree+emdTests = groupTree $$(discover)++prop_iemd_default :: Property+prop_iemd_default = iemdProp defaultEO++prop_orthog_default :: Property+prop_orthog_default = orthogProp defaultEO++edtEO :: EMDOpts Double+edtEO = defaultEO+    { eoSiftCondition = scEnergyDiff 0.01 0.01+                 `SCOr` SCTimes 100+    }++prop_iemd_edt :: Property+prop_iemd_edt = iemdProp edtEO++prop_orthog_edt :: Property+prop_orthog_edt = orthogProp edtEO++sCondEO :: EMDOpts Double+sCondEO = defaultEO+    { eoSiftCondition = SCSCond 10+                 `SCOr` SCTimes 100+    }++prop_iemd_sCond :: Property+prop_iemd_sCond = iemdProp sCondEO++prop_orthog_sCond :: Property+prop_orthog_sCond = orthogProp sCondEO+++iemdProp :: EMDOpts Double -> Property+iemdProp eo = property $ withSize (Range.linear 1 8) $ \(_ :: Proxy n) -> do+    xs <- forAll $ generateData @n+    tripping (CE xs) (emd @_ @_ @(2^n-1) eo . getCE) (Identity . CE . iemd)++orthogProp :: EMDOpts Double -> Property+orthogProp eo = property $ withSize (Range.linear 8 10) $ \(_ :: Proxy n) -> do+    xs   <- forAll $ generateData @n+    let imfs = emdIMFs (emd @_ @_ @(2^n-1) eo xs)+        orthoMatrix =+          [ ((i, j), (x, y), dot x y / sqrt (dot x x * dot y y))+          | (i, x) <- zip indices imfs+          , (j, y) <- zip indices imfs+          , i < j+          ]+        badOrthos = filter (\(_,_,d) -> abs d > 0.5) orthoMatrix+        fracBad :: Double+        fracBad = fromIntegral (length badOrthos)+                / fromIntegral (length orthoMatrix)+    annotateShow orthoMatrix+    annotateShow fracBad+    when (length orthoMatrix < 6) discard+    assert $ fracBad <= 0.5+  where+    indices :: [Int]+    indices = [1..]

@@ -0,0 +1,31 @@+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TemplateHaskell                          #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeInType                               #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}++module Tests.HHT (+    hhtTests+  ) where++import           Data.Functor.Identity+import           Data.Proxy+import           GHC.TypeNats+import           Hedgehog+import           Numeric.EMD+import           Numeric.HHT+import           Test.Tasty+import           Tests.Util+import qualified Hedgehog.Range        as Range++hhtTests :: TestTree+hhtTests = groupTree $$(discover)++prop_ihht :: Property+prop_ihht = property $ withSize (Range.linear 1 5) $ \(_ :: Proxy n) -> do+    xs <- forAll $ generateData @n+    tripping (CE xs) (hht @_ @(2^n-1) defaultEO . getCE) (Identity . CE . ihht)++

@@ -0,0 +1,28 @@+{-# LANGUAGE TypeApplications #-}++module Tests.Spline (+    splineTest+  ) where++import           Data.Maybe+import           Numeric.EMD.Internal.Spline+import           Test.Tasty+import           Test.Tasty.HUnit+import qualified Data.Map                    as M+import qualified Data.Set                    as S++splineTest :: TestTree+splineTest = testCase "Spline Test" $ do+    expected <- map (read @Double) . lines <$> readFile "test-data/sintest.csv"+    roundOut expected @=? roundOut samples+  where+    roundOut :: [Double] -> [Double]+    roundOut = map $ (/ 10e12) . fromInteger . round . (* 10e12)+    spline :: Spline Double+    spline = fromJust+           . makeSpline SENotAKnot+           . M.fromSet sin+           . S.fromList+           $ [0, 1, 2.5, 3.6, 5, 7, 8.1, 10]+    samples :: [Double]+    samples = sampleSpline spline . (/ 4) . fromInteger <$> [0..40]

@@ -0,0 +1,79 @@+{-# LANGUAGE RankNTypes                               #-}+{-# LANGUAGE RecordWildCards                          #-}+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeFamilies                             #-}+{-# LANGUAGE TypeInType                               #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+++module Tests.Util (+    groupTree+  , CloseEnough(..)+  , generateData+  , withSize+  , dot+  ) where++import           Data.Complex+import           Data.Proxy+import           Data.Type.Equality+import           GHC.TypeLits.Compare+import           GHC.TypeNats+import           Hedgehog+import           Hedgehog.Internal.Property+import           Numeric.Natural+import           Statistics.Transform+import           Test.Tasty+import           Test.Tasty.Hedgehog+import qualified Data.Vector.Sized          as V+import qualified Hedgehog.Gen               as Gen+import qualified Hedgehog.Range             as Range++groupTree :: Group -> TestTree+groupTree Group{..} = testGroup (unGroupName groupName)+                                (map (uncurry go) groupProperties)+  where+    go :: PropertyName -> Property -> TestTree+    go n = testProperty (mkName (unPropertyName n))+    mkName = map deUnderscore . drop (length @[] @Char "prop_")+    deUnderscore '_' = ' '+    deUnderscore c   = c++newtype CloseEnough n = CE { getCE :: V.Vector n Double }+  deriving Show++instance KnownNat n => Eq (CloseEnough n) where+    CE x == CE y = ((d `dot` d) / sqrt ((x `dot` x) * (y `dot` y))) < 0.0001+      where+        d = V.zipWith (-) x y++dot :: Num a => V.Vector n a -> V.Vector n a -> a+dot a b = sum $ V.zipWith (*) a b++withSize+    :: Monad m+    => Range Natural+    -> (forall n. (KnownNat n, 1 <= 2^n) => Proxy n -> PropertyT m a)+    -> PropertyT m a+withSize r f = do+    n <- forAll $ Gen.integral r+    case someNatVal n of+      SomeNat (p :: Proxy n) -> do+        LE Refl <- pure $ Proxy @1 %<=? Proxy @(2^n)+        f p++generateData+    :: KnownNat n+    => Gen (V.Vector (2^n) Double)+generateData = fmap (fmap realPart . ifftSized) . V.generateM $ \i ->+    let i' = recip . (+ 1) . fromIntegral $ i+    in  mkPolar <$> Gen.double (Range.exponentialFloat (i' / 10) i')+                <*> Gen.double (Range.constant (-pi) pi)++ifftSized+    :: V.Vector (2^n) (Complex Double)+    -> V.Vector (2^n) (Complex Double)+ifftSized = V.withVectorUnsafe ifft+