emd 0.1.0.0 → 0.1.1.0

raw patch · 5 files changed

+198/−40 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +9/−0
• emd.cabal +3/−2
• src/Numeric/EMD.hs +34/−14
• src/Numeric/EMD/Internal/Spline.hs +23/−24
• src/Numeric/EMD/Unsized.hs +129/−0

@@ -1,6 +1,15 @@ Changelog ========= +Version 0.1.1.0+---------------++*July 25, 2018*++<https://github.com/mstksg/emd/releases/tag/v0.1.1.0>++*   Unsized interface added.+ Version 0.1.0.0 --------------- 

@@ -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 <https://github.com/mstksg/emd#readme> 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

@@ -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

@@ -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 

@@ -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)+