diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,15 @@ Changelog ========= +Version 0.1.1.0 +--------------- + +*July 25, 2018* + + + +* Unsized interface added. + Version 0.1.0.0 --------------- diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 4cb66d27d7d3fcb9f625c87b73448433636262d9476571553b667d1a0f998496 +-- hash: 1007f571ef708a0ed8f93f4737361b9f0bfea8a807a6072e8501352715a67175 name: emd -version: 0.1.0.0 +version: 0.1.1.0 synopsis: Empirical Mode Decomposition (Hilbert-Huang Transform) description: Please see the README on GitHub at category: Math @@ -32,6 +32,7 @@ exposed-modules: Numeric.EMD Numeric.EMD.Internal.Spline + Numeric.EMD.Unsized other-modules: Numeric.EMD.Internal.Tridiagonal Numeric.EMD.Internal.Extrema diff --git a/src/Numeric/EMD.hs b/src/Numeric/EMD.hs --- a/src/Numeric/EMD.hs +++ b/src/Numeric/EMD.hs @@ -24,6 +24,20 @@ -- used to help track down a specific IMF that might be taking more time -- than desired. -- +-- This package uses "sized vectors" as its main interface, to ensure: +-- +-- 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. +-- +-- 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). +-- +-- However, for convenience, "Numeric.EMD.Unsized" is provided with an +-- unsafe unsized interface. +-- module Numeric.EMD ( -- * EMD (Hilbert-Huang Transform) @@ -101,10 +115,16 @@ -- | EMD decomposition (Hilbert-Huang Transform) of a given time series -- with a given sifting stop condition. +-- +-- Takes a sized vector to ensure that: +-- +-- 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) => EMDOpts a - -> SVG.Vector v (n + 2) a - -> EMD v (n + 2) a + -> SVG.Vector v (n + 1) a + -> EMD v (n + 1) a emd eo = runIdentity . emd' (const (pure ())) eo -- | 'emd', but tracing results to stdout as IMFs are found. Useful for @@ -112,8 +132,8 @@ emdTrace :: (VG.Vector v a, KnownNat n, Fractional a, Ord a, MonadIO m) => EMDOpts a - -> SVG.Vector v (n + 2) a - -> m (EMD v (n + 2) a) + -> SVG.Vector v (n + 1) a + -> m (EMD v (n + 1) a) emdTrace = emd' $ \case SRResidual _ -> liftIO $ putStrLn "Residual found." SRIMF _ i -> liftIO $ printf "IMF found (%d iterations)\n" i @@ -121,10 +141,10 @@ -- | 'emd' with a callback for each found IMF. emd' :: (VG.Vector v a, KnownNat n, Fractional a, Ord a, Applicative m) - => (SiftResult v (n + 2) a -> m r) + => (SiftResult v (n + 1) a -> m r) -> EMDOpts a - -> SVG.Vector v (n + 2) a - -> m (EMD v (n + 2) a) + -> SVG.Vector v (n + 1) a + -> m (EMD v (n + 1) a) emd' cb eo = go id where go !imfs !v = cb res *> case res of @@ -136,14 +156,14 @@ -- | 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 iterations + | SRIMF !(SVG.Vector v n a) !Int -- ^ number of 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 + 2) a - -> SiftResult v (n + 2) a + -> SVG.Vector v (n + 1) a + -> SiftResult v (n + 1) a sift EO{..} = go 1 where go !i !v = case sift' eoSplineEnd eoClampEnvelope v of @@ -157,8 +177,8 @@ :: (VG.Vector v a, KnownNat n, Fractional a, Ord a) => SplineEnd -> Bool - -> SVG.Vector v (n + 2) a - -> Maybe (SVG.Vector v (n + 2) a) + -> SVG.Vector v (n + 1) a + -> Maybe (SVG.Vector v (n + 1) a) sift' se cl v = go <$> envelopes se cl v where go (mins, maxs) = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs @@ -170,8 +190,8 @@ :: (VG.Vector v a, KnownNat n, Fractional a, Ord a) => SplineEnd -> Bool - -> SVG.Vector v (n + 2) a - -> Maybe (SVG.Vector v (n + 2) a, SVG.Vector v (n + 2) a) + -> SVG.Vector v (n + 1) a + -> Maybe (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a) envelopes se cl xs = (,) <$> splineAgainst se mins' <*> splineAgainst se maxs' where 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 @@ -43,11 +43,11 @@ | SENatural deriving (Show, Eq, Ord) -data SplineCoef a = SC { scAlpha :: !a -- ^ a - , scBeta :: !a -- ^ b - , scGamma0 :: !a -- ^ y_{i-1} - , scGamma1 :: !a -- ^ y_i - , scDelta :: !a -- ^ x_i - x_{i-1} +data SplineCoef a = SC { _scAlpha :: !a -- ^ a + , _scBeta :: !a -- ^ b + , _scGamma0 :: !a -- ^ y_{i-1} + , _scGamma1 :: !a -- ^ y_i + , _scDelta :: !a -- ^ x_i - x_{i-1} } deriving Show @@ -90,44 +90,43 @@ => SplineEnd -> M.Map a a -- ^ (x, y) -> Maybe (Spline a) -makeSpline se ps = do - (xy0, ps') <- M.minViewWithKey ps - SV.withSizedList (M.toList ps') $ \(xsys :: SV.Vector n (a, a)) -> do +makeSpline se ps = SV.withSizedList (M.toList ps) $ \(xsys :: SV.Vector n (a, a)) -> do Refl <- Proxy @1 `isLE` Proxy @n - let xs, ys :: SV.Vector (n + 1) a - (xs, ys) = SV.unzip $ xy0 `SV.cons` xsys - dxs, dys :: SV.Vector n a + Refl <- Proxy @2 `isLE` Proxy @n + let xs, ys :: SV.Vector n a + (xs, ys) = SV.unzip xsys + dxs, dys :: SV.Vector (n - 1) a dxs = SV.tail xs - SV.init xs - rdxs :: SV.Vector n a + rdxs :: SV.Vector (n - 1) a rdxs = recip dxs - rdxssq :: SV.Vector n a + rdxssq :: SV.Vector (n - 1) a rdxssq = rdxs * rdxs dys = SV.tail ys - SV.init ys dydxssq = dys * rdxssq - mainDiag :: SV.Vector (n - 1) a + mainDiag :: SV.Vector (n - 2) a mainDiag = SV.zipWith (\rdx0 rdx1 -> 2 * ( rdx0 + rdx1 )) (SV.init rdxs) (SV.tail rdxs) - lowerDiag :: SV.Vector (n - 1) a + lowerDiag :: SV.Vector (n - 2) a lowerDiag = SV.take rdxs - upperDiag :: SV.Vector (n - 1) a + upperDiag :: SV.Vector (n - 2) a upperDiag = SV.tail rdxs - rhs :: SV.Vector (n - 1) a + rhs :: SV.Vector (n - 2) a rhs = SV.zipWith (\dydxsq0 dydxsq1 -> 3 * (dydxsq0 + dydxsq1)) (SV.init dydxssq) (SV.tail dydxssq) EE{..} = case se of SENotAKnot -> notAKnot rdxs rdxssq dydxssq SENatural -> natural rdxs dydxssq - solution <- solveTridiagonal (lowerDiag `SV.snoc` eeLower1) - (eeMain0 `SV.cons` mainDiag `SV.snoc` eeMain1) - (eeUpper0 `SV.cons` upperDiag) - (eeRhs0 `SV.cons` rhs `SV.snoc` eeRhs1) - let as :: SV.Vector n a + solution <- solveTridiagonal ( lowerDiag `SV.snoc` eeLower1) + (eeMain0 `SV.cons` mainDiag `SV.snoc` eeMain1 ) + (eeUpper0 `SV.cons` upperDiag ) + (eeRhs0 `SV.cons` rhs `SV.snoc` eeRhs1 ) + let as :: SV.Vector (n - 1) a as = SV.zipWith3 (\k dx dy -> k * dx - dy) (SV.init solution) dxs dys - bs :: SV.Vector n a + bs :: SV.Vector (n - 1) a bs = SV.zipWith3 (\k dx dy -> - k * dx + dy) (SV.tail solution) dxs dys - coefs :: SV.Vector n (a, SplineCoef a) + coefs :: SV.Vector (n - 1) (a, SplineCoef a) coefs = SV.zipWith6 (\x α β γ0 γ1 δ -> (x, SC α β γ0 γ1 δ)) (SV.init xs) as bs (SV.init ys) (SV.tail ys) dxs diff --git a/src/Numeric/EMD/Unsized.hs b/src/Numeric/EMD/Unsized.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/EMD/Unsized.hs @@ -0,0 +1,129 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-} +{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-} + +-- | +-- Module : Numeric.EMD.Unsized +-- Copyright : (c) Justin Le 2018 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- +-- Interface of "Numeric.EMD" re-exported in a non-typesafe "unsized" form. +-- Can be more convenient in certain situations, but "Numeric.EMD" is +-- recommended and preferred. +-- + + +module Numeric.EMD.Unsized ( + -- -- * EMD (Hilbert-Huang Transform) + emd + , emdTrace + , emd' + , EMD(..) + , E.EMDOpts(..), E.defaultEO, E.SiftCondition(..), E.defaultSC, E.SplineEnd(..) + -- -- * Internal + , sift, SiftResult(..) + -- , envelopes + ) where + +import Control.Monad.IO.Class +import Data.Proxy +import Data.Type.Equality +import GHC.TypeLits.Compare +import GHC.TypeNats +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Sized as SVG +import qualified Numeric.EMD as E + +-- | An @'EMD' v a@ is a Hilbert-Huang transform of a time series with +-- items of type @a@ stored in a vector @v@. +data EMD v a = EMD { emdIMFs :: ![v a] + , emdResidual :: !(v a) + } + deriving Show + +-- | The result of a sifting operation. Each sift either yields +-- a residual, or a new IMF. +data SiftResult v a = SRResidual !(v a) + | SRIMF !(v a) !Int -- ^ number of iterations + + +-- | EMD decomposition (Hilbert-Huang Transform) of a given time series +-- with a given sifting stop condition. +-- +-- Returns 'Nothing' if given an empty vector. +-- +-- See 'Numeric.EMD.emd' for a type-safe version with guaruntees on the +-- output vector sizes. +emd :: (VG.Vector v a, Fractional a, Ord a) + => E.EMDOpts a + -> v a + -> Maybe (EMD v a) +emd eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do + Refl <- Proxy @1 `isLE` Proxy @n + pure . convertEMD $ E.emd @_ @_ @(n - 1) eo v' + +-- | 'emd', but tracing results to stdout as IMFs are found. Useful for +-- debugging to see how long each step is taking. +-- +-- Returns 'Nothing' if given an empty vector. +emdTrace + :: (VG.Vector v a, Fractional a, Ord a, MonadIO m) + => E.EMDOpts a + -> v a + -> m (Maybe (EMD v a)) +emdTrace eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> + case Proxy @1 `isLE` Proxy @n of + Nothing -> pure Nothing + Just Refl -> Just . convertEMD <$> E.emdTrace @_ @_ @(n - 1) eo v' + +-- | 'emd' with a callback for each found IMF. +-- +-- Returns 'Nothing' if given an empty vector. +emd' + :: (VG.Vector v a, Ord a, Fractional a, Applicative m) + => (SiftResult v a -> m r) + -> E.EMDOpts a + -> v a + -> m (Maybe (EMD v a)) +emd' cb eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> + case Proxy @1 `isLE` Proxy @n of + Nothing -> pure Nothing + Just Refl -> Just . convertEMD <$> E.emd' @_ @_ @(n - 1) (cb . convertSR) eo v' + +-- emd' cb eo = go id +-- where +-- go !imfs !v = cb res *> case res of +-- SRResidual r -> pure $ EMD (imfs []) r +-- SRIMF v' _ -> go (imfs . (v':)) (v - v') +-- where +-- res = sift eo v + +-- | Iterated sifting process, used to produce either an IMF or a residual. +-- +-- Returns 'Nothing' if given an empty vector. +sift + :: (VG.Vector v a, Fractional a, Ord a) + => E.EMDOpts a + -> v a + -> Maybe (SiftResult v a) +sift eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do + Refl <- Proxy @1 `isLE` Proxy @n + pure $ convertSR . E.sift @_ @_ @(n - 1) eo $ v' + +convertSR :: E.SiftResult v n a -> SiftResult v a +convertSR = \case + E.SRResidual v -> SRResidual $ SVG.fromSized v + E.SRIMF v i -> SRIMF (SVG.fromSized v) i + +convertEMD :: E.EMD v n a -> EMD v a +convertEMD (E.EMD is r) = EMD (SVG.fromSized <$> is) (SVG.fromSized r) +