Unsized.hs

{-# 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)