{-# LANGUAGE NoImplicitPrelude #-}
{- |
Copyright   :  (c) Henning Thielemann 2008
License     :  GPL

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

-}
module Synthesizer.State.Filter.Recursive.MovingAverage
   (sumsStaticInt,
    modulatedFrac,
    ) where

import qualified Synthesizer.State.Signal  as Sig
import qualified Synthesizer.State.Filter.Recursive.Integration as Integration

import qualified Synthesizer.State.Filter.Delay as Delay

import qualified Algebra.Module                as Module
import qualified Algebra.RealField             as RealField

-- import qualified Algebra.Field                 as Field
-- import qualified Algebra.Ring                  as Ring
import qualified Algebra.Additive              as Additive

import NumericPrelude.Base
import NumericPrelude.Numeric



{- |
Like 'Synthesizer.State.Filter.NonRecursive.sums' but in a recursive form.
This needs only linear time (independent of the window size)
but may accumulate rounding errors.

@
ys = xs * (1,0,0,0,-1) \/ (1,-1)
ys * (1,-1) = xs * (1,0,0,0,-1)
ys = xs * (1,0,0,0,-1) + ys * (0,1)
@
-}
{-# INLINE sumsStaticInt #-}
sumsStaticInt :: (Additive.C v) => Int -> Sig.T v -> Sig.T v
sumsStaticInt n xs =
   Integration.run (xs - Delay.staticPos n xs)

{-
staticInt :: (Module.C a v, Additive.C v) => Int -> Sig.T v -> Sig.T v
staticInt n xs =
-}


{-
Sum of a part of a vector with negative sign for reverse order.
It adds from @from@ (inclusively) to @to@ (exclusively),
that is, it sums up @abs (to-from)@ values

sumFromTo :: (Additive.C v) => Int -> Int -> Sig.T v -> v
sumFromTo from to =
   if from <= to
     then          Sig.sum . Sig.take (to-from) . Sig.drop from
     else negate . Sig.sum . Sig.take (from-to) . Sig.drop to
-}

{-# INLINE sumFromToFrac #-}
sumFromToFrac :: (RealField.C a, Module.C a v) => a -> a -> Sig.T v -> v
sumFromToFrac from to xs =
   let (fromInt, fromFrac) = splitFraction from
       (toInt,   toFrac)   = splitFraction to
   in  case compare fromInt toInt of
          EQ -> (to-from) *> Sig.index fromInt xs
          LT ->
            Sig.sum $
            Sig.zipWith id
               (((1-fromFrac) *>) `Sig.cons`
                Sig.replicate (toInt-fromInt-1) id `Sig.append`
                Sig.singleton (toFrac *>)) $
            Sig.drop fromInt xs
          GT ->
            negate $ Sig.sum $
            Sig.zipWith id
               (((1-toFrac) *>) `Sig.cons`
                Sig.replicate (fromInt-toInt-1) id `Sig.append`
                Sig.singleton (fromFrac *>)) $
            Sig.drop toInt xs


{-
            run $
               addNextWeighted (1-toFrac) >>
               replicateM_ (fromInt-toInt-1) addNext >>
               addNextWeighted (fromFrac)

type Accumulator v a =
   WriterT (Dual (Endo v)) (StateT (Sig.T v) Maybe a)

getNext :: Accumulator v a
getNext =
   lift $ StateT $ ListHT.viewL

addAccum :: Additive.C v => v -> Accumulator v ()
addAccum x = tell ((x+) $!)

addNext :: Additive.C v => Accumulator v ()
addNext w =
   addAccum =<< getNext

addNextWeighted :: Module.C a v => a -> Accumulator v ()
addNextWeighted w =
   addAccum . (w *>) =<< getNext
-}

{-
newtype Accumulator v =
   Accumulator ((v, Sig.T v) -> v -> (Sig.T v, v))

addNext :: Additive.C v => Accumulator v
addNext =
   Accumulator $ \(x,xs) s -> (xs, x+s)

addNextWeighted :: Module.C a v => a -> Accumulator v
addNextWeighted a =
   Accumulator $ \(x,xs) s -> (xs, a*>x + s)

bindAccum :: Accumulator v -> Accumulator v -> Accumulator v
bindAccum (Accumulator f) (Accumulator g) =
   Accumulator $ \x s0 ->
      let (ys,s1) = f x s0
      in  maybe s1 () (ListHT.viewL ys)
-}


{- |
Sig.T a must contain only non-negative elements.
-}
{-# INLINE sumDiffsModulated #-}
sumDiffsModulated :: (RealField.C a, Module.C a v) =>
   a -> Sig.T a -> Sig.T v -> Sig.T v
sumDiffsModulated d ds =
   Sig.init .
   -- prevent negative d's since 'drop' cannot restore past values
   Sig.zipWithTails (uncurry sumFromToFrac)
       (Sig.zip (Sig.cons (d+1) ds) (Sig.map (1+) ds)) .
   Sig.cons zero
{-
   Sig.zipWithTails (uncurry sumFromToFrac)
      (Sig.zip (Sig.cons d (Sig.map (subtract 1) ds)) ds)
-}

{-
sumsModulated :: (RealField.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T v -> Sig.T v
sumsModulated maxDInt ds xs =
   let maxD  = fromIntegral maxDInt
       posXs = sumDiffsModulated 0 ds xs
       negXs = sumDiffsModulated maxD (Sig.map (maxD-) ds) (Delay.static maxDInt xs)
   in  Integration.run (posXs - negXs)
-}

{- |
Shift sampling points by a half sample period
in order to preserve signals for window widths below 1.
-}
{-# INLINE sumsModulatedHalf #-}
sumsModulatedHalf :: (RealField.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T v -> Sig.T v
sumsModulatedHalf maxDInt ds xs =
   let maxD  = fromIntegral maxDInt
       d0    = maxD+0.5
       delXs = Delay.staticPos maxDInt xs
       posXs = sumDiffsModulated d0 (Sig.map (d0+) ds) delXs
       negXs = sumDiffsModulated d0 (Sig.map (d0-) ds) delXs
   in  Integration.run (posXs - negXs)

{-# INLINE modulatedFrac #-}
modulatedFrac :: (RealField.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T v -> Sig.T v
modulatedFrac maxDInt ds xs =
   Sig.zipWith (\d y -> recip (2*d) *> y) ds $
   sumsModulatedHalf maxDInt ds xs