diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,25 @@ Changelog ========= +Version 0.1.10.0 +--------------- + +*October 19, 2019* + + + +* 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 --------------- diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -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 diff --git a/src/Numeric/EMD.hs b/src/Numeric/EMD.hs --- a/src/Numeric/EMD.hs +++ b/src/Numeric/EMD.hs @@ -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) diff --git a/src/Numeric/EMD/Internal/Pipe.hs b/src/Numeric/EMD/Internal/Pipe.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/EMD/Internal/Pipe.hs @@ -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 diff --git a/src/Numeric/EMD/Internal/Sift.hs b/src/Numeric/EMD/Internal/Sift.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/EMD/Internal/Sift.hs @@ -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) diff --git a/src/Numeric/EMD/Internal/Spline.hs b/src/Numeric/EMD/Internal/Spline.hs --- a/src/Numeric/EMD/Internal/Spline.hs +++ b/src/Numeric/EMD/Internal/Spline.hs @@ -1,6 +1,6 @@ {-# LANGUAGE ApplicativeDo #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RecordWildCards #-} {-# LANGUAGE ScopedTypeVariables #-} diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -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. diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -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] - diff --git a/test/Tests/EMD.hs b/test/Tests/EMD.hs new file mode 100644 --- /dev/null +++ b/test/Tests/EMD.hs @@ -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..] diff --git a/test/Tests/HHT.hs b/test/Tests/HHT.hs new file mode 100644 --- /dev/null +++ b/test/Tests/HHT.hs @@ -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) + + diff --git a/test/Tests/Spline.hs b/test/Tests/Spline.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Spline.hs @@ -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] diff --git a/test/Tests/Util.hs b/test/Tests/Util.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Util.hs @@ -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 +