{-# LANGUAGE BangPatterns                             #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# OPTIONS_GHC -Wno-orphans                          #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

-- |
-- Module      : Numeric.EMD.Sift
-- Copyright   : (c) Justin Le 2019
-- License     : BSD3
--
-- Maintainer  : [email protected]
-- Stability   : experimental
-- Portability : non-portable
--
-- Tools for creating your own custom sift stopping conditions.
--
-- @since 0.2.0.0
module Numeric.EMD.Sift (
    Sifter(..), SiftResult(..), SingleSift(..), SM
  -- * Sifters
  , defaultSifter
  , siftStdDev
  , siftTimes
  , siftEnergyDiff
  , siftSCond
  , siftAnd
  , siftOr
  -- ** Make Sifters
  , envMean
  , energyDiff
  , normalizeProj
  , siftCauchy
  , siftPairs
  , siftProj
  , siftPairs_
  , siftProj_
  -- * Internal
  , sift, envelopes, rms
  ) where

import           Control.Monad
import           Control.Monad.Trans.Class
import           Control.Monad.Trans.Reader
import           Control.Monad.Trans.State
import           Data.Conduino
import           Data.Conduino.Internal
import           Data.Default.Class
import           Data.Finite
import           Data.Sequence                (Seq(..))
import           GHC.TypeNats
import           Numeric.EMD.Internal
import           Numeric.EMD.Internal.Extrema
import           Numeric.EMD.Internal.Spline
import qualified Data.Conduino.Combinators    as C
import qualified Data.Map                     as M
import qualified Data.Vector.Generic          as VG
import qualified Data.Vector.Generic.Sized    as SVG

-- | @since 0.1.3.0
instance (VG.Vector v a, Fractional a, Ord a) => Default (Sifter v n a) where
    def :: Sifter v n a
def = Sifter v n a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, Fractional a, Ord a) =>
Sifter v n a
defaultSifter

-- | Default 'Sifter'
--
-- @
-- defaultSifter = 'siftStdDev' 0.3 `siftOr` 'siftTimes' 50
-- @
--
-- R package uses @'siftTimes' 20@, Matlab uses no limit
defaultSifter :: (VG.Vector v a, Fractional a, Ord a) => Sifter v n a
defaultSifter :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, Fractional a, Ord a) =>
Sifter v n a
defaultSifter = a -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
(Vector v a, Fractional a, Ord a) =>
a -> Sifter v n a
siftStdDev a
0.3 Sifter v n a -> Sifter v n a -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Sifter v n a -> Sifter v n a -> Sifter v n a
`siftOr` Int -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a. Int -> Sifter v n a
siftTimes Int
50

-- | 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.
siftEnergyDiff
    :: (VG.Vector v a, KnownNat n, Floating a, Ord a)
    => a                -- ^ Threshold for Energy Difference
    -> a                -- ^ Threshold for mean envelope RMS
    -> Sifter v n a
siftEnergyDiff :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a, Ord a) =>
a -> a -> Sifter v n a
siftEnergyDiff a
s a
t = (SingleSift v n a -> SM v n a a) -> a -> Sifter v n a
forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SM v n a b) -> b -> Sifter v n a
siftProj SingleSift v n a -> SM v n a a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, Floating a) =>
SingleSift v n a -> SM v n a a
energyDiff a
s
           Sifter v n a -> Sifter v n a -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Sifter v n a -> Sifter v n a -> Sifter v n a
`siftAnd` (SingleSift v n a -> SM v n a a) -> a -> Sifter v n a
forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SM v n a b) -> b -> Sifter v n a
siftProj SingleSift v n a -> SM v n a a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
SingleSift v n a -> SM v n a a
envMean a
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

-- | Create a sifter that stops after a given fixed number of sifts.
--
-- Useful to use alongside 'siftOr' to set an "upper limit" on the number
-- of sifts.
siftTimes :: Int -> Sifter v n a
siftTimes :: forall (v :: * -> *) (n :: Nat) a. Int -> Sifter v n a
siftTimes Int
n = Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter (Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a)
-> Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall a b. (a -> b) -> a -> b
$ Int -> Pipe (SingleSift v n a) Void Void (SM v n a) ()
forall i o u (m :: * -> *). Int -> Pipe i o u m ()
C.drop (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
forall a b.
Pipe (SingleSift v n a) Void Void (SM v n a) a
-> Pipe (SingleSift v n a) Void Void (SM v n a) b
-> Pipe (SingleSift v n a) Void Void (SM v n a) b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Pipe (SingleSift v n a) Void Void (SM v n a) (SingleSift v n a)
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
forall (f :: * -> *) a. Functor f => f a -> f ()
void Pipe (SingleSift v n a) Void Void (SM v n a) (SingleSift v n a)
forall i o (m :: * -> *). Pipe i o Void m i
awaitSurely

-- | Create a sifter that stops when some projection on 'SingleSift' is
-- smaller than a given threshold.
siftProj
    :: Ord b
    => (SingleSift v n a -> SM v n a b)     -- ^ projection
    -> b                                    -- ^ threshold
    -> Sifter v n a
siftProj :: forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SM v n a b) -> b -> Sifter v n a
siftProj SingleSift v n a -> SM v n a b
p b
t = (SingleSift v n a -> SM v n a Bool) -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
(SingleSift v n a -> SM v n a Bool) -> Sifter v n a
siftProj_ ((SingleSift v n a -> SM v n a Bool) -> Sifter v n a)
-> (SingleSift v n a -> SM v n a Bool) -> Sifter v n a
forall a b. (a -> b) -> a -> b
$ (b -> Bool) -> SM v n a b -> SM v n a Bool
forall a b.
(a -> b)
-> ReaderT (Vector v n a) Identity a
-> ReaderT (Vector v n a) Identity b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (b -> b -> Bool
forall a. Ord a => a -> a -> Bool
<= b
t) (SM v n a b -> SM v n a Bool)
-> (SingleSift v n a -> SM v n a b)
-> SingleSift v n a
-> SM v n a Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SingleSift v n a -> SM v n a b
p

-- | Create a sifter that stops based on some predicate on the initial
-- vector and 'SingleSift' being 'True'.
siftProj_ :: (SingleSift v n a -> SM v n a Bool) -> Sifter v n a
siftProj_ :: forall (v :: * -> *) (n :: Nat) a.
(SingleSift v n a -> SM v n a Bool) -> Sifter v n a
siftProj_ SingleSift v n a -> SM v n a Bool
p = Pipe
  (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
-> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter Pipe
  (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
go
  where
    go :: Pipe
  (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
go = do
      v <- Pipe
  (SingleSift v n a)
  Void
  Void
  (ReaderT (Vector v n a) Identity)
  (SingleSift v n a)
forall i o (m :: * -> *). Pipe i o Void m i
awaitSurely
      r <- lift $ p v
      unless r go

-- | Create a sifter that stops when some projection on two consecutive
-- 'SingleSift's is smaller than a given threshold.
siftPairs
    :: Ord b
    => (SingleSift v n a -> SingleSift v n a -> SM v n a b)
    -> b
    -> Sifter v n a
siftPairs :: forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SingleSift v n a -> SM v n a b)
-> b -> Sifter v n a
siftPairs SingleSift v n a -> SingleSift v n a -> SM v n a b
p b
t = (SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
-> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
(SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
-> Sifter v n a
siftPairs_ ((SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
 -> Sifter v n a)
-> (SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
-> Sifter v n a
forall a b. (a -> b) -> a -> b
$ \SingleSift v n a
x SingleSift v n a
y -> (b -> b -> Bool
forall a. Ord a => a -> a -> Bool
<= b
t) (b -> Bool) -> SM v n a b -> SM v n a Bool
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> SingleSift v n a -> SingleSift v n a -> SM v n a b
p SingleSift v n a
x SingleSift v n a
y

-- | Create a sifter that stops based on some predicate on two consecutive
-- 'SingleSift's being 'True'.
siftPairs_
    :: (SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
    -> Sifter v n a
siftPairs_ :: forall (v :: * -> *) (n :: Nat) a.
(SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
-> Sifter v n a
siftPairs_ SingleSift v n a -> SingleSift v n a -> SM v n a Bool
p = Pipe
  (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
-> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter (Pipe
   (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
 -> Sifter v n a)
-> Pipe
     (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
-> Sifter v n a
forall a b. (a -> b) -> a -> b
$ SingleSift v n a
-> Pipe
     (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
go (SingleSift v n a
 -> Pipe
      (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ())
-> Pipe
     (SingleSift v n a)
     Void
     Void
     (ReaderT (Vector v n a) Identity)
     (SingleSift v n a)
-> Pipe
     (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Pipe
  (SingleSift v n a)
  Void
  Void
  (ReaderT (Vector v n a) Identity)
  (SingleSift v n a)
forall i o (m :: * -> *). Pipe i o Void m i
awaitSurely
  where
    go :: SingleSift v n a
-> Pipe
     (SingleSift v n a) Void Void (ReaderT (Vector v n a) Identity) ()
go SingleSift v n a
s = do
      s' <- Pipe
  (SingleSift v n a)
  Void
  Void
  (ReaderT (Vector v n a) Identity)
  (SingleSift v n a)
forall i o (m :: * -> *). Pipe i o Void m i
awaitSurely
      r  <- lift $ p s s'
      unless r (go s')

-- | Sift based on the "standard deviation test", outlined in original
-- paper.
siftStdDev
    :: forall v n a. (VG.Vector v a, Fractional a, Ord a)
    => a                -- ^ minimal threshold
    -> Sifter v n a
siftStdDev :: forall (v :: * -> *) (n :: Nat) a.
(Vector v a, Fractional a, Ord a) =>
a -> Sifter v n a
siftStdDev = (SingleSift v n a -> SingleSift v n a -> SM v n a a)
-> a -> Sifter v n a
forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SingleSift v n a -> SM v n a b)
-> b -> Sifter v n a
siftPairs ((SingleSift v n a -> SingleSift v n a -> SM v n a a)
 -> a -> Sifter v n a)
-> (SingleSift v n a -> SingleSift v n a -> SM v n a a)
-> a
-> Sifter v n a
forall a b. (a -> b) -> a -> b
$ \(SingleSift Vector v n a
v Vector v n a
_ Vector v n a
_) (SingleSift Vector v n a
v' Vector v n a
_ Vector v n a
_) -> a -> SM v n a a
forall a. a -> ReaderT (Vector v n a) Identity a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (a -> SM v n a a) -> a -> SM v n a a
forall a b. (a -> b) -> a -> b
$
    Vector v n a -> a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, Num a) =>
Vector v n a -> a
SVG.sum ((a -> a -> a) -> Vector v n a -> Vector v n a -> Vector v n a
forall (v :: * -> *) a b c (n :: Nat).
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> Vector v n a -> Vector v n b -> Vector v n c
SVG.zipWith (\a
x a
x' -> (a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
x')a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
xa -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) a -> a -> a
forall a. Num a => a -> a -> a
+ a
eps)) Vector v n a
v Vector v n a
v')
  where
    eps :: a
eps = a
0.0000001

-- | General class of "cauchy-like" sifters: Given a projection function
-- from a 'SingleSift', stop as soon as successive projections become
-- smaller than a given threshold, propertionally.
--
-- Given \(f(x_t)\), stop when:
--
-- \[
--   \frac{(f(x_t) - f(x_{t-1}))^2}{f^2(x_{t-1})} < \delta
-- \]
siftCauchy
    :: (Fractional b, Ord b)
    => (SingleSift v n a -> b)      -- ^ Projection function
    -> b                            -- ^ Threshold \(\delta\)
    -> Sifter v n a
siftCauchy :: forall b (v :: * -> *) (n :: Nat) a.
(Fractional b, Ord b) =>
(SingleSift v n a -> b) -> b -> Sifter v n a
siftCauchy SingleSift v n a -> b
p = (SingleSift v n a -> SingleSift v n a -> SM v n a b)
-> b -> Sifter v n a
forall b (v :: * -> *) (n :: Nat) a.
Ord b =>
(SingleSift v n a -> SingleSift v n a -> SM v n a b)
-> b -> Sifter v n a
siftPairs ((SingleSift v n a -> SingleSift v n a -> SM v n a b)
 -> b -> Sifter v n a)
-> (SingleSift v n a -> SingleSift v n a -> SM v n a b)
-> b
-> Sifter v n a
forall a b. (a -> b) -> a -> b
$ \SingleSift v n a
s SingleSift v n a
s' ->
  let ps :: b
ps  = SingleSift v n a -> b
p SingleSift v n a
s
      ps' :: b
ps' = SingleSift v n a -> b
p SingleSift v n a
s'
      δ :: b
δ   = b
ps' b -> b -> b
forall a. Num a => a -> a -> a
- b
ps
  in  b -> SM v n a b
forall a. a -> ReaderT (Vector v n a) Identity a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (b -> SM v n a b) -> b -> SM v n a b
forall a b. (a -> b) -> a -> b
$ (b
δ b -> b -> b
forall a. Num a => a -> a -> a
* b
δ) b -> b -> b
forall a. Fractional a => a -> a -> a
/ (b
ps b -> b -> b
forall a. Num a => a -> a -> a
* b
ps)

-- | Sift based on the "S-parameter" condition: Stop after a streak @n@ of
-- almost-same numbers of zero crossings and turning points.
siftSCond
    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
    => Int                          -- ^ Streak @n@ to stop on
    -> Sifter v (n + 1) a
siftSCond :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
Int -> Sifter v (n + 1) a
siftSCond Int
n = Pipe (SingleSift v (n + 1) a) Void Void (SM v (n + 1) a) ()
-> Sifter v (n + 1) a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter (Pipe (SingleSift v (n + 1) a) Void Void (SM v (n + 1) a) ()
 -> Sifter v (n + 1) a)
-> Pipe (SingleSift v (n + 1) a) Void Void (SM v (n + 1) a) ()
-> Sifter v (n + 1) a
forall a b. (a -> b) -> a -> b
$ (SingleSift v (n + 1) a -> Int)
-> Pipe (SingleSift v (n + 1) a) Int Void (SM v (n + 1) a) Void
forall i o u (m :: * -> *). (i -> o) -> Pipe i o u m u
C.map (Vector v (n + 1) a -> Int
forall {a} {v :: * -> *} {n :: Nat}.
(Vector v a, KnownNat n, Fractional a, Ord a) =>
Vector v (n + 1) a -> Int
crossCount (Vector v (n + 1) a -> Int)
-> (SingleSift v (n + 1) a -> Vector v (n + 1) a)
-> SingleSift v (n + 1) a
-> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SingleSift v (n + 1) a -> Vector v (n + 1) a
forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssResult)
                    Pipe (SingleSift v (n + 1) a) Int Void (SM v (n + 1) a) Void
-> Pipe Int Void Void (SM v (n + 1) a) ()
-> Pipe (SingleSift v (n + 1) a) Void Void (SM v (n + 1) a) ()
forall (m :: * -> *) a b u v c r.
Monad m =>
Pipe a b u m v -> Pipe b c v m r -> Pipe a c u m r
.| Int -> Pipe Int (Seq Int) Void (SM v (n + 1) a) Void
forall i u (m :: * -> *). Int -> Pipe i (Seq i) u m u
C.consecutive Int
n
                    Pipe Int (Seq Int) Void (SM v (n + 1) a) Void
-> Pipe (Seq Int) Void Void (SM v (n + 1) a) ()
-> Pipe Int Void Void (SM v (n + 1) a) ()
forall (m :: * -> *) a b u v c r.
Monad m =>
Pipe a b u m v -> Pipe b c v m r -> Pipe a c u m r
.| (Seq Int -> Maybe (Seq Int, Int))
-> Pipe (Seq Int) (Seq Int, Int) Void (SM v (n + 1) a) Void
forall (t :: * -> *) i o u (m :: * -> *).
Foldable t =>
(i -> t o) -> Pipe i o u m u
C.concatMap Seq Int -> Maybe (Seq Int, Int)
pick
                    Pipe (Seq Int) (Seq Int, Int) Void (SM v (n + 1) a) Void
-> Pipe (Seq Int, Int) Void Void (SM v (n + 1) a) ()
-> Pipe (Seq Int) Void Void (SM v (n + 1) a) ()
forall (m :: * -> *) a b u v c r.
Monad m =>
Pipe a b u m v -> Pipe b c v m r -> Pipe a c u m r
.| ((Seq Int, Int) -> Bool)
-> Pipe (Seq Int, Int) Void Void (SM v (n + 1) a) ()
forall i o u (m :: * -> *). (i -> Bool) -> Pipe i o u m ()
C.dropWhile (Seq Int, Int) -> Bool
forall {t :: * -> *} {b}.
(Foldable t, Num b, Ord b) =>
(t b, b) -> Bool
notGood
  where
    pick :: Seq Int -> Maybe (Seq Int, Int)
pick Seq Int
Empty      = Maybe (Seq Int, Int)
forall a. Maybe a
Nothing
    pick (Seq Int
xs :|> Int
x) = (Seq Int
xs, Int
x) (Seq Int, Int) -> Maybe () -> Maybe (Seq Int, Int)
forall a b. a -> Maybe b -> Maybe a
forall (f :: * -> *) a b. Functor f => a -> f b -> f a
<$ Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Seq Int -> Int
forall a. Seq a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Seq Int
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
    notGood :: (t b, b) -> Bool
notGood (t b
xs, b
x) = (b -> Bool) -> t b -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((b -> b -> Bool
forall a. Ord a => a -> a -> Bool
<= b
1) (b -> Bool) -> (b -> b) -> b -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b
forall a. Num a => a -> a
abs (b -> b) -> (b -> b) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b -> b
forall a. Num a => a -> a -> a
subtract b
x) t b
xs
    crossCount :: Vector v (n + 1) a -> Int
crossCount Vector v (n + 1) a
xs = Map (Finite (n + 1)) a -> Int
forall k a. Map k a -> Int
M.size Map (Finite (n + 1)) a
mins Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Map (Finite (n + 1)) a -> Int
forall k a. Map k a -> Int
M.size Map (Finite (n + 1)) a
maxs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
crosses
      where
        (Map (Finite (n + 1)) a
mins, Map (Finite (n + 1)) a
maxs) = Vector v (n + 1) a
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
forall (v :: * -> *) (n :: Nat) a.
(Vector v a, KnownNat n, Fractional a, Ord a) =>
Vector v (n + 1) a
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
extrema Vector v (n + 1) a
xs
        crosses :: Int
crosses = (Int, Maybe Bool) -> Int
forall a b. (a, b) -> a
fst ((Int, Maybe Bool) -> Int)
-> ((a -> StateT (Int, Maybe Bool) Identity ())
    -> (Int, Maybe Bool))
-> (a -> StateT (Int, Maybe Bool) Identity ())
-> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (StateT (Int, Maybe Bool) Identity ()
 -> (Int, Maybe Bool) -> (Int, Maybe Bool))
-> (Int, Maybe Bool)
-> StateT (Int, Maybe Bool) Identity ()
-> (Int, Maybe Bool)
forall a b c. (a -> b -> c) -> b -> a -> c
flip StateT (Int, Maybe Bool) Identity ()
-> (Int, Maybe Bool) -> (Int, Maybe Bool)
forall s a. State s a -> s -> s
execState (Int
0, Maybe Bool
forall a. Maybe a
Nothing) (StateT (Int, Maybe Bool) Identity () -> (Int, Maybe Bool))
-> ((a -> StateT (Int, Maybe Bool) Identity ())
    -> StateT (Int, Maybe Bool) Identity ())
-> (a -> StateT (Int, Maybe Bool) Identity ())
-> (Int, Maybe Bool)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((a -> StateT (Int, Maybe Bool) Identity ())
 -> Vector v (n + 1) a -> StateT (Int, Maybe Bool) Identity ())
-> Vector v (n + 1) a
-> (a -> StateT (Int, Maybe Bool) Identity ())
-> StateT (Int, Maybe Bool) Identity ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip (a -> StateT (Int, Maybe Bool) Identity ())
-> Vector v (n + 1) a -> StateT (Int, Maybe Bool) Identity ()
forall (m :: * -> *) (v :: * -> *) a b (n :: Nat).
(Monad m, Vector v a) =>
(a -> m b) -> Vector v n a -> m ()
SVG.mapM_ Vector v (n + 1) a
xs ((a -> StateT (Int, Maybe Bool) Identity ()) -> Int)
-> (a -> StateT (Int, Maybe Bool) Identity ()) -> Int
forall a b. (a -> b) -> a -> b
$ \a
x -> ((Int, Maybe Bool) -> (Int, Maybe Bool))
-> StateT (Int, Maybe Bool) Identity ()
forall (m :: * -> *) s. Monad m => (s -> s) -> StateT s m ()
modify (((Int, Maybe Bool) -> (Int, Maybe Bool))
 -> StateT (Int, Maybe Bool) Identity ())
-> ((Int, Maybe Bool) -> (Int, Maybe Bool))
-> StateT (Int, Maybe Bool) Identity ()
forall a b. (a -> b) -> a -> b
$ \(!Int
i, !Maybe Bool
y) ->
          let xPos :: Bool
xPos = a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0
              i' :: Int
i'   = case Maybe Bool
y of
                       Maybe Bool
Nothing -> Int
i
                       Just Bool
y'
                         | Bool
xPos Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== Bool
y' -> Int
i
                         | Bool
otherwise  -> Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
          in  (Int
i', Bool -> Maybe Bool
forall a. a -> Maybe a
Just Bool
xPos)

-- | Combine two sifters in "or" fashion: The final sifter will complete
-- when /either/ sifter completes.
siftOr :: Sifter v n a -> Sifter v n a -> Sifter v n a
siftOr :: forall (v :: * -> *) (n :: Nat) a.
Sifter v n a -> Sifter v n a -> Sifter v n a
siftOr (Sifter Pipe (SingleSift v n a) Void Void (SM v n a) ()
p) (Sifter Pipe (SingleSift v n a) Void Void (SM v n a) ()
q) = Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter (Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a)
-> Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall a b. (a -> b) -> a -> b
$ Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
forall (m :: * -> *) i u a.
Monad m =>
Pipe i Void u m a -> Pipe i Void u m a -> Pipe i Void u m a
altSink Pipe (SingleSift v n a) Void Void (SM v n a) ()
p Pipe (SingleSift v n a) Void Void (SM v n a) ()
q
infixr 2 `siftOr`

-- | Combine two sifters in "and" fashion: The final sifter will complete
-- when /both/ sifters complete.
siftAnd :: Sifter v n a -> Sifter v n a -> Sifter v n a
siftAnd :: forall (v :: * -> *) (n :: Nat) a.
Sifter v n a -> Sifter v n a -> Sifter v n a
siftAnd (Sifter Pipe (SingleSift v n a) Void Void (SM v n a) ()
p) (Sifter Pipe (SingleSift v n a) Void Void (SM v n a) ()
q) = Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall (v :: * -> *) (n :: Nat) a.
Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
Sifter (Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a)
-> Pipe (SingleSift v n a) Void Void (SM v n a) () -> Sifter v n a
forall a b. (a -> b) -> a -> b
$ Pipe (SingleSift v n a) Void Void (SM v n a) (() -> ())
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
forall (m :: * -> *) i u a b.
Monad m =>
Pipe i Void u m (a -> b) -> Pipe i Void u m a -> Pipe i Void u m b
zipSink (() -> ()
forall a. a -> a
id (() -> ())
-> Pipe (SingleSift v n a) Void Void (SM v n a) ()
-> Pipe (SingleSift v n a) Void Void (SM v n a) (() -> ())
forall a b.
a
-> Pipe (SingleSift v n a) Void Void (SM v n a) b
-> Pipe (SingleSift v n a) Void Void (SM v n a) a
forall (f :: * -> *) a b. Functor f => a -> f b -> f a
<$ Pipe (SingleSift v n a) Void Void (SM v n a) ()
p) Pipe (SingleSift v n a) Void Void (SM v n a) ()
q
infixr 3 `siftAnd`

-- | Project the root mean square of the mean of the maximum and minimum
-- envelopes.
envMean
    :: (VG.Vector v a, KnownNat n, Floating a)
    => SingleSift v n a
    -> SM v n a a
envMean :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
SingleSift v n a -> SM v n a a
envMean SingleSift{Vector v n a
ssResult :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssResult :: Vector v n a
ssMinEnv :: Vector v n a
ssMaxEnv :: Vector v n a
ssMaxEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssMinEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
..} = a -> ReaderT (Vector v n a) Identity a
forall a. a -> ReaderT (Vector v n a) Identity a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (a -> ReaderT (Vector v n a) Identity a)
-> a -> ReaderT (Vector v n a) Identity a
forall a b. (a -> b) -> a -> b
$
    Vector v n a -> a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
Vector v n a -> a
rms (Vector v n a -> a) -> Vector v n a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> Vector v n a -> Vector v n a -> Vector v n a
forall (v :: * -> *) a b c (n :: Nat).
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> Vector v n a -> Vector v n b -> Vector v n c
SVG.zipWith (\a
x a
y -> (a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
y) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2) Vector v n a
ssMinEnv Vector v n a
ssMaxEnv

-- | Project the /square root/ of the "Energy difference".
energyDiff
    :: (VG.Vector v a, Floating a)
    => SingleSift v n a
    -> SM v n a a
energyDiff :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, Floating a) =>
SingleSift v n a -> SM v n a a
energyDiff SingleSift{Vector v n a
ssResult :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssMaxEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssMinEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssResult :: Vector v n a
ssMinEnv :: Vector v n a
ssMaxEnv :: Vector v n a
..} = do
    v0 <- ReaderT (Vector v n a) Identity (Vector v n a)
forall (m :: * -> *) r. Monad m => ReaderT r m r
ask
    pure . sqrt . abs . SVG.sum
         $ SVG.zipWith (\a
x a
c -> a
c a -> a -> a
forall a. Num a => a -> a -> a
* (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
c)) v0 ssResult

-- | Given a "projection function" (like 'envMean' or 'energyDiff'),
-- re-scale the result based on the RMS of the original signal.
normalizeProj
    :: (VG.Vector v a, KnownNat n, Floating a)
    => (SingleSift v n a -> SM v n a a)
    -> (SingleSift v n a -> SM v n a a)
normalizeProj :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
(SingleSift v n a -> SM v n a a) -> SingleSift v n a -> SM v n a a
normalizeProj SingleSift v n a -> SM v n a a
f SingleSift v n a
ss = do
    v0 <- (Vector v n a -> a) -> SM v n a a
forall (m :: * -> *) r a. Monad m => (r -> a) -> ReaderT r m a
asks Vector v n a -> a
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
Vector v n a -> a
rms
    r  <- f ss
    pure $ r / v0

-- | Get the root mean square of a vector
rms :: (VG.Vector v a, KnownNat n, Floating a) => SVG.Vector v n a -> a
rms :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Floating a) =>
Vector v n a -> a
rms Vector v n a
xs = a -> a
forall a. Floating a => a -> a
sqrt (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> a -> Vector v n a -> a
forall (v :: * -> *) b a (n :: Nat).
Vector v b =>
(a -> b -> a) -> a -> Vector v n b -> a
SVG.foldl' (\a
s a
x -> a
s a -> a -> a
forall a. Num a => a -> a -> a
+ a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
x) a
0 Vector v n a
xs a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Vector v n a -> Int
forall (v :: * -> *) (n :: Nat) a.
KnownNat n =>
Vector v n a -> Int
SVG.length Vector v n a
xs)


-- | 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 v (n + 1) a
    -> SVG.Vector v (n + 1) a
    -> SiftResult v (n + 1) a
sift :: forall (v :: * -> *) (n :: Nat) a.
(Vector v a, KnownNat n, Floating a, Ord a) =>
EMDOpts v (n + 1) a -> Vector v (n + 1) a -> SiftResult v (n + 1) a
sift EO{Maybe BoundaryHandler
SplineEnd a
Sifter v (n + 1) a
eoSifter :: Sifter v (n + 1) a
eoSplineEnd :: SplineEnd a
eoBoundaryHandler :: Maybe BoundaryHandler
eoBoundaryHandler :: forall (v :: * -> *) (n :: Nat) a.
EMDOpts v n a -> Maybe BoundaryHandler
eoSifter :: forall (v :: * -> *) (n :: Nat) a. EMDOpts v n a -> Sifter v n a
eoSplineEnd :: forall (v :: * -> *) (n :: Nat) a. EMDOpts v n a -> SplineEnd a
..} Vector v (n + 1) a
v0 = case StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) ()
-> (Int, Vector v (n + 1) a)
-> Either (Vector v (n + 1) a) (Int, Vector v (n + 1) a)
forall (m :: * -> *) s a. Monad m => StateT s m a -> s -> m s
execStateT (Pipe
  ()
  Void
  (ZonkAny 0)
  (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
  ()
-> StateT
     (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) ()
forall (m :: * -> *) u a. Monad m => Pipe () Void u m a -> m a
runPipe Pipe
  ()
  Void
  (ZonkAny 0)
  (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
  ()
sifterPipe) (Int
0, Vector v (n + 1) a
v0) of
    Left  Vector v (n + 1) a
v        -> Vector v (n + 1) a -> SiftResult v (n + 1) a
forall (v :: * -> *) (n :: Nat) a. Vector v n a -> SiftResult v n a
SRResidual Vector v (n + 1) a
v
    Right (!Int
i, !Vector v (n + 1) a
v) -> Vector v (n + 1) a -> Int -> SiftResult v (n + 1) a
forall (v :: * -> *) (n :: Nat) a.
Vector v n a -> Int -> SiftResult v n a
SRIMF Vector v (n + 1) a
v Int
i
  where
    sifterPipe :: Pipe
  ()
  Void
  (ZonkAny 0)
  (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
  ()
sifterPipe = StateT
  (Int, Vector v (n + 1) a)
  (Either (Vector v (n + 1) a))
  (SingleSift v (n + 1) a)
-> Pipe
     ()
     (SingleSift v (n + 1) a)
     (ZonkAny 0)
     (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
     Void
forall (m :: * -> *) o i u a. Monad m => m o -> Pipe i o u m a
C.repeatM StateT
  (Int, Vector v (n + 1) a)
  (Either (Vector v (n + 1) a))
  (SingleSift v (n + 1) a)
go
              Pipe
  ()
  (SingleSift v (n + 1) a)
  (ZonkAny 0)
  (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
  Void
-> Pipe
     (SingleSift v (n + 1) a)
     Void
     Void
     (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
     ()
-> Pipe
     ()
     Void
     (ZonkAny 0)
     (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
     ()
forall (m :: * -> *) a b u v c r.
Monad m =>
Pipe a b u m v -> Pipe b c v m r -> Pipe a c u m r
.| (forall x.
 ReaderT (Vector v (n + 1) a) Identity x
 -> StateT
      (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) x)
-> Pipe
     (SingleSift v (n + 1) a)
     Void
     Void
     (ReaderT (Vector v (n + 1) a) Identity)
     ()
-> Pipe
     (SingleSift v (n + 1) a)
     Void
     Void
     (StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)))
     ()
forall (m :: * -> *) (n :: * -> *) i o u a.
(Monad m, Monad n) =>
(forall x. m x -> n x) -> Pipe i o u m a -> Pipe i o u n a
hoistPipe
                    (x
-> StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) x
forall a.
a
-> StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (x
 -> StateT
      (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) x)
-> (ReaderT (Vector v (n + 1) a) Identity x -> x)
-> ReaderT (Vector v (n + 1) a) Identity x
-> StateT (Int, Vector v (n + 1) a) (Either (Vector v (n + 1) a)) x
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ReaderT (Vector v (n + 1) a) Identity x -> Vector v (n + 1) a -> x
forall r a. Reader r a -> r -> a
`runReader` Vector v (n + 1) a
v0))
                    (Sifter v (n + 1) a
-> Pipe
     (SingleSift v (n + 1) a)
     Void
     Void
     (ReaderT (Vector v (n + 1) a) Identity)
     ()
forall (v :: * -> *) (n :: Nat) a.
Sifter v n a -> Pipe (SingleSift v n a) Void Void (SM v n a) ()
sPipe Sifter v (n + 1) a
eoSifter)
    go :: StateT
  (Int, Vector v (n + 1) a)
  (Either (Vector v (n + 1) a))
  (SingleSift v (n + 1) a)
go = ((Int, Vector v (n + 1) a)
 -> Either
      (Vector v (n + 1) a)
      (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a)))
-> StateT
     (Int, Vector v (n + 1) a)
     (Either (Vector v (n + 1) a))
     (SingleSift v (n + 1) a)
forall s (m :: * -> *) a. (s -> m (a, s)) -> StateT s m a
StateT (((Int, Vector v (n + 1) a)
  -> Either
       (Vector v (n + 1) a)
       (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a)))
 -> StateT
      (Int, Vector v (n + 1) a)
      (Either (Vector v (n + 1) a))
      (SingleSift v (n + 1) a))
-> ((Int, Vector v (n + 1) a)
    -> Either
         (Vector v (n + 1) a)
         (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a)))
-> StateT
     (Int, Vector v (n + 1) a)
     (Either (Vector v (n + 1) a))
     (SingleSift v (n + 1) a)
forall a b. (a -> b) -> a -> b
$ \(!Int
i, !Vector v (n + 1) a
v) ->
      case SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (SingleSift v (n + 1) a)
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (SingleSift v (n + 1) a)
sift' SplineEnd a
eoSplineEnd Maybe BoundaryHandler
eoBoundaryHandler Vector v (n + 1) a
v of
        Maybe (SingleSift v (n + 1) a)
Nothing                -> Vector v (n + 1) a
-> Either
     (Vector v (n + 1) a)
     (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a))
forall a b. a -> Either a b
Left Vector v (n + 1) a
v
        Just ss :: SingleSift v (n + 1) a
ss@SingleSift{Vector v (n + 1) a
ssResult :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssMaxEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssMinEnv :: forall (v :: * -> *) (n :: Nat) a. SingleSift v n a -> Vector v n a
ssResult :: Vector v (n + 1) a
ssMinEnv :: Vector v (n + 1) a
ssMaxEnv :: Vector v (n + 1) a
..} -> (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a))
-> Either
     (Vector v (n + 1) a)
     (SingleSift v (n + 1) a, (Int, Vector v (n + 1) a))
forall a b. b -> Either a b
Right (SingleSift v (n + 1) a
ss, (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Vector v (n + 1) a
ssResult))

-- | 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' :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (SingleSift v (n + 1) a)
sift' SplineEnd a
se Maybe BoundaryHandler
bh Vector v (n + 1) a
v = do
    (mins, maxs) <- SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (Vector v (n + 1) a, Vector v (n + 1) a)
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (Vector v (n + 1) a, Vector v (n + 1) a)
envelopes SplineEnd a
se Maybe BoundaryHandler
bh Vector v (n + 1) a
v
    pure SingleSift
      { ssResult = SVG.zipWith3 (\a
x a
mi a
ma -> a
x a -> a -> a
forall a. Num a => a -> a -> a
- (a
mi a -> a -> a
forall a. Num a => a -> a -> a
+ a
ma)a -> a -> a
forall a. Fractional a => a -> a -> a
/a
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 :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Maybe BoundaryHandler
-> Vector v (n + 1) a
-> Maybe (Vector v (n + 1) a, Vector v (n + 1) a)
envelopes SplineEnd a
se Maybe BoundaryHandler
bh Vector v (n + 1) a
xs = do
    Bool -> Maybe () -> Maybe ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Maybe BoundaryHandler
bh Maybe BoundaryHandler -> Maybe BoundaryHandler -> Bool
forall a. Eq a => a -> a -> Bool
== BoundaryHandler -> Maybe BoundaryHandler
forall a. a -> Maybe a
Just BoundaryHandler
BHClamp) (Maybe () -> Maybe ()) -> Maybe () -> Maybe ()
forall a b. (a -> b) -> a -> b
$ do
      Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Map (Finite (n + 1)) a -> Int
forall k a. Map k a -> Int
M.size Map (Finite (n + 1)) a
mins Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1)
      Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Map (Finite (n + 1)) a -> Int
forall k a. Map k a -> Int
M.size Map (Finite (n + 1)) a
maxs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1)
    (,) (Vector v (n + 1) a
 -> Vector v (n + 1) a -> (Vector v (n + 1) a, Vector v (n + 1) a))
-> Maybe (Vector v (n + 1) a)
-> Maybe
     (Vector v (n + 1) a -> (Vector v (n + 1) a, Vector v (n + 1) a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> SplineEnd a
-> Map Int a
-> Map (Finite (n + 1)) a
-> Maybe (Vector v (n + 1) a)
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Map Int a -> Map (Finite n) a -> Maybe (Vector v n a)
splineAgainst SplineEnd a
se Map Int a
emin Map (Finite (n + 1)) a
mins
        Maybe
  (Vector v (n + 1) a -> (Vector v (n + 1) a, Vector v (n + 1) a))
-> Maybe (Vector v (n + 1) a)
-> Maybe (Vector v (n + 1) a, Vector v (n + 1) a)
forall a b. Maybe (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> SplineEnd a
-> Map Int a
-> Map (Finite (n + 1)) a
-> Maybe (Vector v (n + 1) a)
forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Map Int a -> Map (Finite n) a -> Maybe (Vector v n a)
splineAgainst SplineEnd a
se Map Int a
emax Map (Finite (n + 1)) a
maxs
  where
    -- minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]
    (Map (Finite (n + 1)) a
mins,Map (Finite (n + 1)) a
maxs) = Vector v (n + 1) a
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
forall (v :: * -> *) (n :: Nat) a.
(Vector v a, KnownNat n, Fractional a, Ord a) =>
Vector v (n + 1) a
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
extrema Vector v (n + 1) a
xs
    (Map Int a
emin,Map Int a
emax) = case Maybe BoundaryHandler
bh of
      Maybe BoundaryHandler
Nothing  -> (Map Int a, Map Int a)
forall a. Monoid a => a
mempty
      Just BoundaryHandler
bh' -> Vector v (n + 1) a
-> BoundaryHandler
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
-> (Map Int a, Map Int a)
forall (v :: * -> *) (n :: Nat) a.
(Vector v a, KnownNat n) =>
Vector v (n + 1) a
-> BoundaryHandler
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
-> (Map Int a, Map Int a)
extendExtrema Vector v (n + 1) a
xs BoundaryHandler
bh' (Map (Finite (n + 1)) a
mins,Map (Finite (n + 1)) a
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 :: forall (v :: * -> *) (n :: Nat) a.
(Vector v a, KnownNat n) =>
Vector v (n + 1) a
-> BoundaryHandler
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
-> (Map Int a, Map Int a)
extendExtrema Vector v (n + 1) a
xs = \case
    BoundaryHandler
BHClamp     -> (Map Int a, Map Int a)
-> (Map (Finite (n + 1)) a, Map (Finite (n + 1)) a)
-> (Map Int a, Map Int a)
forall a b. a -> b -> a
const (Map Int a
firstLast, Map Int a
firstLast)
    BoundaryHandler
BHSymmetric -> \(Map (Finite (n + 1)) a
mins, Map (Finite (n + 1)) a
maxs) ->
      let addFirst :: (Map Int a, Map Int a)
addFirst = case (Maybe (Finite (n + 1), Map Int a)
flippedMin, Maybe (Finite (n + 1), Map Int a)
flippedMax) of
              (Maybe (Finite (n + 1), Map Int a)
Nothing      , Maybe (Finite (n + 1), Map Int a)
Nothing      ) -> (Map Int a, Map Int a)
forall a. Monoid a => a
mempty
              -- first point is local maximum
              (Just (Finite (n + 1)
_,Map Int a
mn)  , Maybe (Finite (n + 1), Map Int a)
Nothing      ) -> (Map Int a
mn        , Map Int a
firstPoint)
              -- first point is local minimum
              (Maybe (Finite (n + 1), Map Int a)
Nothing      , Just (Finite (n + 1)
_,Map Int a
mx)  ) -> (Map Int a
firstPoint, Map Int a
mx        )
              (Just (Finite (n + 1)
mni,Map Int a
mn), Just (Finite (n + 1)
mxi,Map Int a
mx))
                | Finite (n + 1)
mni Finite (n + 1) -> Finite (n + 1) -> Bool
forall a. Ord a => a -> a -> Bool
< Finite (n + 1)
mxi                  -> (Map Int a
mn        , Map Int a
firstPoint)
                | Bool
otherwise                  -> (Map Int a
firstPoint, Map Int a
mx        )
            where
              flippedMin :: Maybe (Finite (n + 1), Map Int a)
flippedMin = (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a)
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Map (Finite (n + 1)) a -> Maybe (Finite (n + 1), a)
forall k a. Map k a -> Maybe (k, a)
M.lookupMin Map (Finite (n + 1)) a
mins) (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), Map Int a))
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> a -> b
$ \(Finite (n + 1)
minIx, a
minVal) ->
                (Finite (n + 1)
minIx, Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton (Int -> Int
forall a. Num a => a -> a
negate (Finite (n + 1) -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Finite (n + 1)
minIx)) a
minVal)
              flippedMax :: Maybe (Finite (n + 1), Map Int a)
flippedMax = (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a)
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Map (Finite (n + 1)) a -> Maybe (Finite (n + 1), a)
forall k a. Map k a -> Maybe (k, a)
M.lookupMin Map (Finite (n + 1)) a
maxs) (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), Map Int a))
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> a -> b
$ \(Finite (n + 1)
maxIx, a
maxVal) ->
                (Finite (n + 1)
maxIx, Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton (Int -> Int
forall a. Num a => a -> a
negate (Finite (n + 1) -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Finite (n + 1)
maxIx)) a
maxVal)
          addLast :: (Map Int a, Map Int a)
addLast = case (Maybe (Finite (n + 1), Map Int a)
flippedMin, Maybe (Finite (n + 1), Map Int a)
flippedMax) of
              (Maybe (Finite (n + 1), Map Int a)
Nothing      , Maybe (Finite (n + 1), Map Int a)
Nothing      ) -> (Map Int a, Map Int a)
forall a. Monoid a => a
mempty
              -- last point is local maximum
              (Just (Finite (n + 1)
_,Map Int a
mn)  , Maybe (Finite (n + 1), Map Int a)
Nothing      ) -> (Map Int a
mn        , Map Int a
lastPoint )
              -- last point is local minimum
              (Maybe (Finite (n + 1), Map Int a)
Nothing      , Just (Finite (n + 1)
_,Map Int a
mx)  ) -> (Map Int a
lastPoint , Map Int a
mx        )
              (Just (Finite (n + 1)
mni,Map Int a
mn), Just (Finite (n + 1)
mxi,Map Int a
mx))
                | Finite (n + 1)
mni Finite (n + 1) -> Finite (n + 1) -> Bool
forall a. Ord a => a -> a -> Bool
> Finite (n + 1)
mxi                  -> (Map Int a
mn        , Map Int a
lastPoint )
                | Bool
otherwise                  -> (Map Int a
lastPoint , Map Int a
mx        )
            where
              flippedMin :: Maybe (Finite (n + 1), Map Int a)
flippedMin = (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a)
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Map (Finite (n + 1)) a -> Maybe (Finite (n + 1), a)
forall k a. Map k a -> Maybe (k, a)
M.lookupMax Map (Finite (n + 1)) a
mins) (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), Map Int a))
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> a -> b
$ \(Finite (n + 1)
minIx, a
minVal) ->
                (Finite (n + 1)
minIx, Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton (Int -> Int
extendSym (Finite (n + 1) -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Finite (n + 1)
minIx)) a
minVal)
              flippedMax :: Maybe (Finite (n + 1), Map Int a)
flippedMax = (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a)
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), a) -> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Map (Finite (n + 1)) a -> Maybe (Finite (n + 1), a)
forall k a. Map k a -> Maybe (k, a)
M.lookupMax Map (Finite (n + 1)) a
maxs) (((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
 -> Maybe (Finite (n + 1), Map Int a))
-> ((Finite (n + 1), a) -> (Finite (n + 1), Map Int a))
-> Maybe (Finite (n + 1), Map Int a)
forall a b. (a -> b) -> a -> b
$ \(Finite (n + 1)
maxIx, a
maxVal) ->
                (Finite (n + 1)
maxIx, Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton (Int -> Int
extendSym (Finite (n + 1) -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Finite (n + 1)
maxIx)) a
maxVal)
      in  (Map Int a, Map Int a)
addFirst (Map Int a, Map Int a)
-> (Map Int a, Map Int a) -> (Map Int a, Map Int a)
forall a. Monoid a => a -> a -> a
`mappend` (Map Int a, Map Int a)
addLast
  where
    lastIx :: Int
lastIx = Finite n -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Finite n -> Int) -> Finite n -> Int
forall a b. (a -> b) -> a -> b
$ forall a. Bounded a => a
maxBound @(Finite n)
    firstPoint :: Map Int a
firstPoint = Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton Int
0 (Vector v (1 + n) a -> a
forall (v :: * -> *) (n :: Nat) a.
Vector v a =>
Vector v (1 + n) a -> a
SVG.head Vector v (n + 1) a
Vector v (1 + n) a
xs)
    lastPoint :: Map Int a
lastPoint  = Int -> a -> Map Int a
forall k a. k -> a -> Map k a
M.singleton Int
lastIx (Vector v (n + 1) a -> a
forall (v :: * -> *) (n :: Nat) a.
Vector v a =>
Vector v (n + 1) a -> a
SVG.last Vector v (n + 1) a
xs)
    firstLast :: Map Int a
firstLast  = Map Int a
firstPoint Map Int a -> Map Int a -> Map Int a
forall a. Monoid a => a -> a -> a
`mappend` Map Int a
lastPoint
    extendSym :: Int -> Int
extendSym Int
i = Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
lastIx Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
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 :: forall (v :: * -> *) a (n :: Nat).
(Vector v a, KnownNat n, Fractional a, Ord a) =>
SplineEnd a
-> Map Int a -> Map (Finite n) a -> Maybe (Vector v n a)
splineAgainst SplineEnd a
se Map Int a
ext = (Spline a -> Vector v n a)
-> Maybe (Spline a) -> Maybe (Vector v n a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Spline a -> Vector v n a
forall {n :: Nat} {v :: * -> *} {a}.
(KnownNat n, Vector v a, Fractional a, Ord a) =>
Spline a -> Vector v n a
go
                     (Maybe (Spline a) -> Maybe (Vector v n a))
-> (Map (Finite n) a -> Maybe (Spline a))
-> Map (Finite n) a
-> Maybe (Vector v n a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SplineEnd a -> Map a a -> Maybe (Spline a)
forall a.
(Ord a, Fractional a) =>
SplineEnd a -> Map a a -> Maybe (Spline a)
makeSpline SplineEnd a
se
                     (Map a a -> Maybe (Spline a))
-> (Map (Finite n) a -> Map a a)
-> Map (Finite n) a
-> Maybe (Spline a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map a a -> Map a a -> Map a a
forall a. Monoid a => a -> a -> a
mappend ((Int -> a) -> Map Int a -> Map a a
forall k1 k2 a. (k1 -> k2) -> Map k1 a -> Map k2 a
M.mapKeysMonotonic Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Map Int a
ext)
                     (Map a a -> Map a a)
-> (Map (Finite n) a -> Map a a) -> Map (Finite n) a -> Map a a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Finite n -> a) -> Map (Finite n) a -> Map a a
forall k1 k2 a. (k1 -> k2) -> Map k1 a -> Map k2 a
M.mapKeysMonotonic Finite n -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral
  where
    go :: Spline a -> Vector v n a
go Spline a
spline = (Finite n -> a) -> Vector v n a
forall (v :: * -> *) (n :: Nat) a.
(KnownNat n, Vector v a) =>
(Finite n -> a) -> Vector v n a
SVG.generate (Spline a -> a -> a
forall a. (Fractional a, Ord a) => Spline a -> a -> a
sampleSpline Spline a
spline (a -> a) -> (Finite n -> a) -> Finite n -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Finite n -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral)