diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,16 @@ Changelog ========= +Version 0.1.5.0 +--------------- + +*August 31, 2019* + + + +* Add `NFData` instance for `EMD`, `HHTLine`, and `HTT` +* Add `iemd`, inverting `emd`. + Version 0.1.4.0 --------------- diff --git a/bench/bench.hs b/bench/bench.hs new file mode 100644 --- /dev/null +++ b/bench/bench.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeInType #-} +{-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-} + +import Control.DeepSeq +import Control.Exception +import Criterion.Main +import Data.Char +import Data.Complex +import Data.Maybe +import GHC.TypeNats +import Numeric.EMD +import Text.Printf +import qualified Data.Vector as UV +import qualified Data.Vector.Sized as V +import qualified Numeric.FFT as FFT +import qualified System.Random.MWC as MWC + +main :: IO () +main = do + g <- MWC.initialize + . UV.fromList + . map (fromIntegral . ord) + $ "hello world" + + test256 <- evaluate . force =<< generateData @8 g + test1024 <- evaluate . force =<< generateData @10 g + test4096 <- evaluate . force =<< generateData @12 g + test16384 <- evaluate . force =<< generateData @14 g + + itest256 <- evaluate . force $ emd defaultEO test256 + itest1024 <- evaluate . force $ emd defaultEO test1024 + itest4096 <- evaluate . force $ emd defaultEO test4096 + itest16384 <- evaluate . force $ emd defaultEO test16384 + + let imfs256 = length . emdIMFs $ itest256 + imfs1024 = length . emdIMFs $ itest1024 + imfs4096 = length . emdIMFs $ itest4096 + imfs16384 = length . emdIMFs $ itest16384 + + defaultMainWith defaultConfig [ + bgroup "emd" + [ bench (printf "256 (%d imfs)" imfs256 ) $ nf (emd defaultEO) test256 + , bench (printf "1024 (%d imfs)" imfs1024 ) $ nf (emd defaultEO) test1024 + , bench (printf "4096 (%d imfs)" imfs4096 ) $ nf (emd defaultEO) test4096 + , bench (printf "16384 (%d imfs)" imfs16384) $ nf (emd defaultEO) test16384 + ] + , bgroup "iemd" + [ bench "256" $ nf iemd itest256 + , bench "1024" $ nf iemd itest1024 + , bench "4096" $ nf iemd itest4096 + , bench "16384" $ nf iemd itest16384 + ] + ] + +generateData + :: KnownNat n + => MWC.GenIO + -> IO (V.Vector (2^n) Double) +generateData g = fmap (fmap realPart . ifftSized) . V.generateM $ \i -> + let i' = recip . (+ 1) . fromIntegral $ i + in (:+) <$> MWC.uniformR (-i', i') g + <*> MWC.uniformR (-i', i') g + where + + +ifftSized + :: KnownNat n + => V.Vector (2^n) (Complex Double) + -> V.Vector (2^n) (Complex Double) +ifftSized = fromJust + . V.fromList + . FFT.ifft + . V.toList diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -1,11 +1,13 @@ --- This file has been generated from package.yaml by hpack version 0.28.2. +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.31.2. -- -- see: https://github.com/sol/hpack -- --- hash: cae340ee4756e6a3c4e72c728779352322ea841f2e7230892d25b139b413e5f1 +-- hash: 8f9970f61b5024a25a75d52714da895f6afaacbafe209e2021934b996a00d91a name: emd -version: 0.1.4.0 +version: 0.1.5.0 synopsis: Empirical Mode Decomposition and Hilbert-Huang Transform description: Please see the README on GitHub at category: Math @@ -18,10 +20,9 @@ license-file: LICENSE tested-with: GHC >= 8.2 build-type: Simple -cabal-version: >= 1.10 extra-source-files: - CHANGELOG.md README.md + CHANGELOG.md test-data/sintest.csv source-repository head @@ -44,6 +45,7 @@ , binary , containers , data-default-class + , deepseq , finite-typelits , ghc-typelits-knownnat , ghc-typelits-natnormalise @@ -66,4 +68,24 @@ , base >=4.10 && <5 , containers , emd + default-language: Haskell2010 + +benchmark emd-bench + type: exitcode-stdio-1.0 + main-is: bench.hs + other-modules: + Paths_emd + hs-source-dirs: + bench + ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N -O2 + build-depends: + base >=4.10 && <5 + , criterion + , deepseq + , emd + , ghc-typelits-knownnat + , mwc-random + , pure-fft + , vector + , vector-sized default-language: Haskell2010 diff --git a/src/Numeric/EMD.hs b/src/Numeric/EMD.hs --- a/src/Numeric/EMD.hs +++ b/src/Numeric/EMD.hs @@ -43,6 +43,7 @@ emd , emdTrace , emd' + , iemd , EMD(..) , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..) -- * Internal @@ -50,11 +51,13 @@ , 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.TypeNats import Numeric.EMD.Internal.Extrema @@ -157,12 +160,15 @@ -- with @n@ items of type @a@ stored in a vector @v@. -- -- The component-wise sum of 'emdIMFs' and 'emdResidual' should yield --- exactly the original series. +-- exactly the original series (see 'iemd'). data EMD v n a = EMD { emdIMFs :: ![SVG.Vector v n a] , emdResidual :: !(SVG.Vector v n a) } deriving (Show, Generic, Eq, Ord) +-- | @since 0.1.5.0 +instance NFData (v a) => NFData (EMD v n a) + -- | @since 0.1.3.0 instance (VG.Vector v a, KnownNat n, Bi.Binary (v a)) => Bi.Binary (EMD v n a) where put EMD{..} = Bi.put (SVG.fromSized <$> emdIMFs) @@ -211,6 +217,17 @@ SRIMF v' _ -> go (imfs . (v':)) (v - v') where res = sift eo v + +-- | Collapse an 'EMD' back into its original time series. Should be +-- a left-inverse to 'emd': using 'iemd' on the result of 'emd' should give +-- back the original vector. +-- +-- @since 0.1.5.0 +iemd + :: (VG.Vector v a, Num a) + => 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. diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -24,6 +24,12 @@ -- 'hhtEmd'. See "Numeric.EMD" for information on why this module uses -- "sized vectors", and how to convert unsized vectors to sized vectors. -- +-- Note that the Hilbert Transform implementation in this module is +-- slightly naive and is essentially O(n^2) on the length of the vector. +-- However, computation time for the full Hilbert-Huang Transform is +-- typically dominated by Empirical Mode Docomposition, which is +-- approximately O(n). +-- -- @since 0.1.2.0 module Numeric.HHT ( @@ -44,6 +50,7 @@ , hilbertMagFreq ) where +import Control.DeepSeq import Data.Complex import Data.Finite import Data.Fixed @@ -79,6 +86,9 @@ Just hlFreqs <- SVG.toSized <$> Bi.get pure HHTLine{..} +-- | @since 0.1.5.0 +instance NFData (v a) => NFData (HHTLine v n a) + -- | A Hilbert-Huang Transform. An @'HHT' v n a@ is a Hilbert-Huang -- transform of an @n@-item time series of items of type @a@ represented -- using vector @v@. @@ -90,6 +100,9 @@ -- | @since 0.1.3.0 instance (VG.Vector v a, KnownNat n, Bi.Binary (v a)) => Bi.Binary (HHT v n a) +-- | @since 0.1.5.0 +instance NFData (v a) => NFData (HHT v n a) + -- | Directly compute the Hilbert-Huang transform of a given time series. -- Essentially is a composition of 'hhtEmd' and 'emd'. See 'hhtEmd' for -- a more flexible version. @@ -210,7 +223,7 @@ -- | Returns the dominant frequency (frequency with largest magnitude -- contribution) at each time step. --- +-- -- @since 0.1.4.0 dominantFreq :: forall v n a. (VG.Vector v a, KnownNat n, Ord a)