emd 0.1.10.0 → 0.2.0.0

raw patch · 9 files changed

+495/−413 lines, 9 filesdep +conduinoPVP ok

version bump matches the API change (PVP)

Dependencies added: conduino

API changes (from Hackage documentation)

Files

• CHANGELOG.md +10/−0
• emd.cabal +6/−3
• src/Numeric/EMD.hs +36/−12
• src/Numeric/EMD/Internal.hs +61/−0
• src/Numeric/EMD/Internal/Sift.hs +0/−384
• src/Numeric/EMD/Internal/Spline.hs +1/−1
• src/Numeric/EMD/Sift.hs +365/−0
• src/Numeric/HHT.hs +3/−3
• test/Tests/EMD.hs +13/−10

@@ -1,6 +1,16 @@ Changelog ========= +Version 0.2.0.0+---------------++*October 30, 2019*++<https://github.com/mstksg/emd/releases/tag/v0.2.0.0>++*   Sift condition system totally revamped, allowing for custom sift+    conditions.+ Version 0.1.10.0 --------------- 

@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 5230a7607e7d3f26121daa88b982f8f58087a236adc506e45eadf1a038581d7a+-- hash: 3facd98354d5dd814f0bea76e951d633090222a995f013f939679dc36a60f586  name:           emd-version:        0.1.10.0+version:        0.2.0.0 synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform description:    Empirical Mode decomposition and Hilbert-Huang Transform in pure                 Haskell.@@ -34,11 +34,12 @@   exposed-modules:       Numeric.EMD       Numeric.EMD.Internal.Spline+      Numeric.EMD.Sift       Numeric.HHT   other-modules:+      Numeric.EMD.Internal       Numeric.EMD.Internal.Extrema       Numeric.EMD.Internal.Pipe-      Numeric.EMD.Internal.Sift       Numeric.EMD.Internal.Tridiagonal       Numeric.HHT.Internal.FFT   hs-source-dirs:@@ -49,6 +50,7 @@     , base >=4.11 && <5     , binary     , carray+    , conduino     , containers     , data-default-class     , deepseq@@ -88,6 +90,7 @@     , tasty-hedgehog     , tasty-hunit     , typelits-witnesses+    , vector     , vector-sized   default-language: Haskell2010 

@@ -8,12 +8,13 @@ {-# LANGUAGE TypeApplications                         #-} {-# LANGUAGE TypeInType                               #-} {-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -Wno-orphans                          #-} {-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-} {-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}  -- | -- Module      : Numeric.EMD--- Copyright   : (c) Justin Le 2018+-- Copyright   : (c) Justin Le 2019 -- License     : BSD3 -- -- Maintainer  : justin@jle.im@@ -38,7 +39,6 @@ -- '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 (   -- * Empirical Mode Decomposition     emd@@ -49,8 +49,8 @@   -- ** Configuration   , EMDOpts(..), defaultEO   , BoundaryHandler(..)-  , SiftCondition(..), SiftProjection(..), defaultSC-  , scEnergyDiff+  , Sifter+  , defaultSifter   , SplineEnd(..)   -- * Internal   , sift, SiftResult(..)@@ -59,16 +59,40 @@  import           Control.DeepSeq import           Control.Monad.IO.Class+import           Data.Default.Class import           Data.Functor.Identity import           Data.List-import           GHC.Generics              (Generic)+import           GHC.Generics                (Generic) import           GHC.TypeNats-import           Numeric.EMD.Internal.Sift+import           Numeric.EMD.Internal+import           Numeric.EMD.Internal.Spline+import           Numeric.EMD.Sift import           Text.Printf-import qualified Data.Binary               as Bi-import qualified Data.Vector.Generic       as VG-import qualified Data.Vector.Generic.Sized as SVG+import qualified Data.Binary                 as Bi+import qualified Data.Vector.Generic         as VG+import qualified Data.Vector.Generic.Sized   as SVG +-- | Default 'EMDOpts'+--+-- Note: If you immediately use this and set 'eoSifter', then @v@ will be+-- ambiguous.  Explicitly set @v@ with type applications to appease GHC+--+-- @+-- 'defaultEO' @(Data.Vector.Vector)+--    { eoSifter = scTimes 100+--    }+-- @+defaultEO :: (VG.Vector v a, Fractional a, Ord a) => EMDOpts v n a+defaultEO = EO { eoSifter          = defaultSifter+               , eoSplineEnd       = SENatural+               , eoBoundaryHandler = Just BHSymmetric+               }++-- | @since 0.1.3.0+instance (VG.Vector v a, Fractional a, Ord a) => Default (EMDOpts v n a) where+    def = defaultEO++ -- | 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@. --@@ -100,7 +124,7 @@ --     the input vector -- 2.  We provide a vector of size of at least one. emd :: (VG.Vector v a, KnownNat n, Floating a, Ord a)-    => EMDOpts a+    => EMDOpts v (n + 1) a     -> SVG.Vector v (n + 1) a     -> EMD v (n + 1) a emd eo = runIdentity . emd' (const (pure ())) eo@@ -109,7 +133,7 @@ -- debugging to see how long each step is taking. emdTrace     :: (VG.Vector v a, KnownNat n, Floating a, Ord a, MonadIO m)-    => EMDOpts a+    => EMDOpts v (n + 1) a     -> SVG.Vector v (n + 1) a     -> m (EMD v (n + 1) a) emdTrace = emd' $ \case@@ -120,7 +144,7 @@ emd'     :: (VG.Vector v a, KnownNat n, Floating a, Ord a, Applicative m)     => (SiftResult v (n + 1) a -> m r)-    -> EMDOpts a+    -> EMDOpts v (n + 1) a     -> SVG.Vector v (n + 1) a     -> m (EMD v (n + 1) a) emd' cb eo = go id

@@ -0,0 +1,61 @@+{-# LANGUAGE DeriveGeneric #-}++module Numeric.EMD.Internal (+    EMDOpts(..)+  , BoundaryHandler(..)+  , Sifter(..)+  , SM+  , SingleSift(..)+  ) where++import           Control.Monad.Trans.Reader+import           Data.Conduino+import           Data.Void+import           GHC.Generics+import           Numeric.EMD.Internal.Spline+import qualified Data.Binary                 as Bi+import qualified Data.Vector.Generic.Sized   as SVG++-- | Options for EMD composition.+data EMDOpts v n a = EO+    { eoSifter          :: Sifter v n a           -- ^ stop condition for sifting+    , eoSplineEnd       :: SplineEnd a            -- ^ end conditions for envelope splines+    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary+    }+  deriving (Generic)++-- | 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++-- | Result of a single sift+data SingleSift v n a = SingleSift+    { ssResult :: !(SVG.Vector v n a)+    , ssMinEnv :: !(SVG.Vector v n a)+    , ssMaxEnv :: !(SVG.Vector v n a)+    }++-- | Monad where 'Sifter' actions live.  The reader parameter is the+-- "original vector".+type SM v n a = Reader (SVG.Vector v n a)++-- | A sift stopping condition.+--+-- It is a 'Pipe' consumer that takes single sift step results upstream and+-- terminates with '()' as soon as it is satisfied with the latest sift+-- step.+--+-- Use combinators like 'siftOr' and 'siftAnd' to combine sifters, and the+-- various sifters in "Numeric.EMD.Sift" to create sifters from commonly+-- established ones or new ones from scratch.+--+-- @since 0.2.0.0+newtype Sifter v n a = Sifter { sPipe :: Pipe (SingleSift v n a) Void Void (SM v n a) () }+

@@ -1,384 +0,0 @@-{-# 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)

@@ -12,7 +12,7 @@  -- | -- Module      : Numeric.EMD.Internal.Spline--- Copyright   : (c) Justin Le 2018+-- Copyright   : (c) Justin Le 2019 -- License     : BSD3 -- -- Maintainer  : justin@jle.im

@@ -0,0 +1,365 @@+{-# LANGUAGE BangPatterns                             #-}+{-# LANGUAGE GADTs                                    #-}+{-# LANGUAGE LambdaCase                               #-}+{-# LANGUAGE RankNTypes                               #-}+{-# LANGUAGE RecordWildCards                          #-}+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeInType                               #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -Wno-orphans                          #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}++-- |+-- Module      : Numeric.EMD.Sift+-- Copyright   : (c) Justin Le 2019+-- License     : BSD3+--+-- Maintainer  : justin@jle.im+-- Stability   : experimental+-- Portability : non-portable+--+-- Tools for creating your own custom sift stopping conditions.+--+-- @since 0.2.0.0+module Numeric.EMD.Sift (+    Sifter(..), SiftResult(..), SingleSift(..), SM+  -- * Sifters+  , defaultSifter+  , siftStdDev+  , siftTimes+  , siftEnergyDiff+  , siftSCond+  , siftAnd+  , siftOr+  -- ** Make Sifters+  , envMean+  , energyDiff+  , normalizeProj+  , siftCauchy+  , siftPairs+  , siftProj+  , siftPairs_+  , siftProj_+  -- * Internal+  , sift, envelopes, rms+  ) where++import           Control.Monad+import           Control.Monad.Trans.Class+import           Control.Monad.Trans.Reader+import           Control.Monad.Trans.State+import           Data.Conduino+import           Data.Conduino.Internal+import           Data.Default.Class+import           Data.Finite+import           Data.Sequence                (Seq(..))+import           GHC.TypeNats+import           Numeric.EMD.Internal+import           Numeric.EMD.Internal.Extrema+import           Numeric.EMD.Internal.Spline+import qualified Data.Conduino.Combinators    as C+import qualified Data.Map                     as M+import qualified Data.Vector.Generic          as VG+import qualified Data.Vector.Generic.Sized    as SVG++-- | @since 0.1.3.0+instance (VG.Vector v a, Fractional a, Ord a) => Default (Sifter v n a) where+    def = defaultSifter++-- | Default 'Sifter'+--+-- @+-- defaultSifter = 'siftStdDev' 0.3 `siftOr` 'siftTimes' 50+-- @+--+-- R package uses @'siftTimes' 20@, Matlab uses no limit+defaultSifter :: (VG.Vector v a, Fractional a, Ord a) => Sifter v n a+defaultSifter = siftStdDev 0.3 `siftOr` siftTimes 50++-- | 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.+siftEnergyDiff+    :: (VG.Vector v a, KnownNat n, Floating a, Ord a)+    => a                -- ^ Threshold for Energy Difference+    -> a                -- ^ Threshold for mean envelope RMS+    -> Sifter v n a+siftEnergyDiff s t = siftProj energyDiff s+           `siftAnd` siftProj envMean 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++-- | Create a sifter that stops after a given fixed number of sifts.+--+-- Useful to use alongside 'siftOr' to set an "upper limit" on the number+-- of sifts.+siftTimes :: Int -> Sifter v n a+siftTimes n = Sifter $ C.drop (n - 1) >> void awaitSurely++-- | Create a sifter that stops when some projection on 'SingleSift' is+-- smaller than a given threshold.+siftProj+    :: Ord b+    => (SingleSift v n a -> SM v n a b)     -- ^ projection+    -> b                                    -- ^ threshold+    -> Sifter v n a+siftProj p t = siftProj_ $ fmap (<= t) . p++-- | Create a sifter that stops based on some predicate on the initial+-- vector and 'SingleSift' being 'True'.+siftProj_ :: (SingleSift v n a -> SM v n a Bool) -> Sifter v n a+siftProj_ p = Sifter go+  where+    go = do+      v <- awaitSurely+      r <- lift $ p v+      unless r go++-- | Create a sifter that stops when some projection on two consecutive+-- 'SingleSift's is smaller than a given threshold.+siftPairs+    :: Ord b+    => (SingleSift v n a -> SingleSift v n a -> SM v n a b)+    -> b+    -> Sifter v n a+siftPairs p t = siftPairs_ $ \x y -> (<= t) <$> p x y++-- | Create a sifter that stops based on some predicate on two consecutive+-- 'SingleSift's being 'True'.+siftPairs_+    :: (SingleSift v n a -> SingleSift v n a -> SM v n a Bool)+    -> Sifter v n a+siftPairs_ p = Sifter $ go =<< awaitSurely+  where+    go s = do+      s' <- awaitSurely+      r  <- lift $ p s s'+      unless r (go s')++-- | Sift based on the "standard deviation test", outlined in original+-- paper.+siftStdDev+    :: forall v n a. (VG.Vector v a, Fractional a, Ord a)+    => a                -- ^ minimal threshold+    -> Sifter v n a+siftStdDev = siftPairs $ \(SingleSift v _ _) (SingleSift v' _ _) -> pure $+    SVG.sum (SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v')+  where+    eps = 0.0000001++-- | General class of "cauchy-like" sifters: Given a projection function+-- from a 'SingleSift', stop as soon as successive projections become+-- smaller than a given threshold, propertionally.+--+-- Given \(f(x_t)\), stop when:+--+-- \[+--   \frac{(f(x_t) - f(x_{t-1}))^2}{f^2(x_{t-1})} < \delta+-- \]+siftCauchy+    :: (Fractional b, Ord b)+    => (SingleSift v n a -> b)      -- ^ Projection function+    -> b                            -- ^ Threshold \(\delta\)+    -> Sifter v n a+siftCauchy p = siftPairs $ \s s' ->+  let ps  = p s+      ps' = p s'+      δ   = ps' - ps+  in  pure $ (δ * δ) / (ps * ps)++-- | Sift based on the "S-parameter" condition: Stop after a streak @n@ of+-- almost-same numbers of zero crossings and turning points.+siftSCond+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)+    => Int                          -- ^ Streak @n@ to stop on+    -> Sifter v (n + 1) a+siftSCond n = Sifter $ C.map (crossCount . ssResult)+                    .| C.consecutive n+                    .| C.concatMap pick+                    .| C.dropWhile notGood+  where+    pick Empty      = Nothing+    pick (xs :|> x) = (xs, x) <$ guard (length xs == (n - 1))+    notGood (xs, x) = all ((<= 1) . abs . subtract x) xs+    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)++-- | Combine two sifters in "or" fashion: The final sifter will complete+-- when /either/ sifter completes.+siftOr :: Sifter v n a -> Sifter v n a -> Sifter v n a+siftOr (Sifter p) (Sifter q) = Sifter $ altSink p q+infixr 2 `siftOr`++-- | Combine two sifters in "and" fashion: The final sifter will complete+-- when /both/ sifters complete.+siftAnd :: Sifter v n a -> Sifter v n a -> Sifter v n a+siftAnd (Sifter p) (Sifter q) = Sifter $ zipSink (id <$ p) q+infixr 3 `siftAnd`++-- | Project the root mean square of the mean of the maximum and minimum+-- envelopes.+envMean+    :: (VG.Vector v a, KnownNat n, Floating a)+    => SingleSift v n a+    -> SM v n a a+envMean SingleSift{..} = pure $+    rms $ SVG.zipWith (\x y -> (x + y) / 2) ssMinEnv ssMaxEnv++-- | Project the /square root/ of the "Energy difference".+energyDiff+    :: (VG.Vector v a, Floating a)+    => SingleSift v n a+    -> SM v n a a+energyDiff SingleSift{..} = do+    v0 <- ask+    pure . sqrt . abs . SVG.sum+         $ SVG.zipWith (\x c -> c * (x - c)) v0 ssResult++-- | Given a "projection function" (like 'envMean' or 'energyDiff'),+-- re-scale the result based on the RMS of the original signal.+normalizeProj+    :: (VG.Vector v a, KnownNat n, Floating a)+    => (SingleSift v n a -> SM v n a a)+    -> (SingleSift v n a -> SM v n a a)+normalizeProj f ss = do+    v0 <- asks rms+    r  <- f ss+    pure $ r / v0++-- | Get the root mean square of a vector+rms :: (VG.Vector v a, KnownNat n, Floating a) => SVG.Vector v n a -> a+rms xs = sqrt $ SVG.foldl' (\s x -> s + x*x) 0 xs / fromIntegral (SVG.length xs)+++-- | 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 v (n + 1) 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 = C.repeatM go+              .| hoistPipe+                    (pure . (`runReader` v0))+                    (sPipe eoSifter)+    go = StateT $ \(!i, !v) ->+      case sift' eoSplineEnd eoBoundaryHandler v of+        Nothing                -> Left v+        Just ss@SingleSift{..} -> Right (ss, (i + 1, ssResult))++-- | 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+      { ssResult = 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)

@@ -14,7 +14,7 @@  -- | -- Module      : Numeric.HHT--- Copyright   : (c) Justin Le 2018+-- Copyright   : (c) Justin Le 2019 -- License     : BSD3 -- -- Maintainer  : justin@jle.im@@ -43,7 +43,7 @@   , expectedFreq, dominantFreq   , foldFreq   -- ** Options-  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)+  , EMDOpts(..), defaultEO, BoundaryHandler(..), defaultSifter, SplineEnd(..)   -- * Hilbert transforms (internal usage)   , hilbert   , hilbertIm@@ -125,7 +125,7 @@ -- Essentially is a composition of 'hhtEmd' and 'emd'.  See 'hhtEmd' for -- a more flexible version. hht :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)-    => EMDOpts a+    => EMDOpts v (n + 1) a     -> SVG.Vector v (n + 1) a     -> HHT v n a hht eo = hhtEmd . emd eo

@@ -1,3 +1,4 @@+{-# LANGUAGE RankNTypes                               #-} {-# LANGUAGE ScopedTypeVariables                      #-} {-# LANGUAGE TemplateHaskell                          #-} {-# LANGUAGE TypeApplications                         #-}@@ -16,8 +17,10 @@ import           GHC.TypeNats import           Hedgehog import           Numeric.EMD+import           Numeric.EMD.Sift import           Test.Tasty import           Tests.Util+import qualified Data.Vector           as UV import qualified Hedgehog.Range        as Range  emdTests :: TestTree@@ -29,10 +32,10 @@ prop_orthog_default :: Property prop_orthog_default = orthogProp defaultEO -edtEO :: EMDOpts Double-edtEO = defaultEO-    { eoSiftCondition = scEnergyDiff 0.01 0.01-                 `SCOr` SCTimes 100+edtEO :: KnownNat n => EMDOpts UV.Vector n Double+edtEO = (defaultEO @UV.Vector)+    { eoSifter = siftEnergyDiff 0.01 0.01+        `siftOr` siftTimes 100     }  prop_iemd_edt :: Property@@ -41,10 +44,10 @@ prop_orthog_edt :: Property prop_orthog_edt = orthogProp edtEO -sCondEO :: EMDOpts Double-sCondEO = defaultEO-    { eoSiftCondition = SCSCond 10-                 `SCOr` SCTimes 100+sCondEO :: KnownNat n => EMDOpts UV.Vector (n + 1) Double+sCondEO = (defaultEO @UV.Vector)+    { eoSifter = siftSCond 10+        `siftOr` siftTimes 100     }  prop_iemd_sCond :: Property@@ -54,12 +57,12 @@ prop_orthog_sCond = orthogProp sCondEO  -iemdProp :: EMDOpts Double -> Property+iemdProp :: (forall n. KnownNat n => EMDOpts UV.Vector (n + 1) 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 :: (forall n. KnownNat n => EMDOpts UV.Vector (n + 1) 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)