{-# OPTIONS -fno-implicit-prelude #-}
{- |
ToDo:
use AffineSpace instead of Module for the particular interpolation types,
since affine combinations assert reconstruction of constant functions.
They are more natural for interpolation of internal control parameters.
However, how can cubic interpolation expressed by affine combinations
without divisions?
-}
module Synthesizer.Generic.Interpolation where

import qualified Synthesizer.Generic.Control      as Ctrl
import qualified Synthesizer.Generic.SampledValue as Sample
import qualified Synthesizer.Generic.Signal       as SigG

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 Algebra.Additive(zero)
import Algebra.Module((*>))
import Data.Maybe (fromMaybe)

import Control.Monad.State (StateT(StateT), evalStateT, ap, )
import Control.Applicative (Applicative(pure, (<*>)), (<$>), liftA2, )
import Synthesizer.ApplicativeUtility (liftA4, )
import Synthesizer.Utility (affineComb, )

import PreludeBase
import NumericPrelude


{- | interpolation as needed for resampling -}
data T sig t y =
  Cons {
    number :: Int,  -- interpolation requires a total number of 'number'
    offset :: Int,  -- interpolation requires 'offset' values before the current
    func   :: t -> sig y -> y
  }


{-* Interpolation with various padding methods -}

{-# INLINE zeroPad #-}
zeroPad :: (RealField.C t, Sample.C y, SigG.C sig) =>
   (T sig t y -> t -> sig y -> a) ->
   y -> T sig t y -> t -> sig y -> a
zeroPad interpolate z ip phase x =
   let (phInt, phFrac) = splitFraction phase
   in  interpolate ip phFrac
          (delayPad z (offset ip - phInt) (SigG.append x (SigG.repeat z)))

{-# INLINE constantPad #-}
constantPad :: (RealField.C t, Sample.C y, SigG.C sig) =>
   (T sig t y -> t -> sig y -> a) ->
   T sig t y -> t -> sig y -> a
constantPad interpolate ip phase x =
   let (phInt, phFrac) = splitFraction phase
       xPad =
          do (xFirst,_) <- SigG.viewL x
             return (delayPad xFirst (offset ip - phInt) (SigG.extendConstant x))
   in  interpolate ip phFrac
          (fromMaybe SigG.empty xPad)


{- |
Only for finite input signals.
-}
{-# INLINE cyclicPad #-}
cyclicPad :: (RealField.C t, Sample.C y, SigG.C sig) =>
   (T sig t y -> t -> sig y -> a) ->
   T sig t y -> t -> sig y -> a
cyclicPad interpolate ip phase x =
   let (phInt, phFrac) = splitFraction phase
   in  interpolate ip phFrac
          (SigG.drop (mod (phInt - offset ip) (SigG.length x)) (SigG.cycle x))

{- |
The extrapolation may miss some of the first and some of the last points
-}
{-# INLINE extrapolationPad #-}
extrapolationPad :: (RealField.C t, Sample.C y, SigG.C sig) =>
   (T sig t y -> t -> sig y -> a) ->
   T sig t y -> t -> sig y -> a
extrapolationPad interpolate ip phase =
   interpolate ip (phase - fromIntegral (offset ip))
{-
  This example shows pikes, although there shouldn't be any:
   plotList (take 100 $ interpolate (Zero (0::Double)) ipCubic (-0.9::Double) (repeat 0.03) [1,0,1,0.8])
-}


{-* Interpolation of multiple values with various padding methods -}

{-# INLINE skip #-}
skip :: (RealField.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> (t, sig y) -> (t, sig y)
skip ip (phase0, x0) =
   let (n, frac) = splitFraction phase0
       (m, x1) = SigG.dropMarginRem (number ip) n x0
   in  (fromIntegral m + frac, x1)

{-# INLINE single #-}
single :: (RealField.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> t -> sig y -> y
single ip phase0 x0 =
   uncurry (func ip) $ skip ip (phase0, x0)
--   curry (uncurry (func ip) . skip ip)
{-
GNUPlot.plotFunc [] (GNUPlot.linearScale 1000 (0,2)) (\t -> single linear (t::Double) [0,4,1::Double])
-}


{-* Interpolation of multiple values with various padding methods -}

{- | All values of frequency control must be non-negative. -}
{-# INLINE multiRelative #-}
multiRelative ::
   (RealField.C t, Sample.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> t -> sig y -> sig t -> sig y
multiRelative ip phase0 x0 =
   SigG.crochetL
      (\freq pos ->
          let (phase,x) = skip ip pos
          in  Just (func ip phase x, (phase+freq,x)))
      (phase0,x0)


{-# INLINE multiRelativeZeroPad #-}
multiRelativeZeroPad ::
   (RealField.C t, Sample.C t, Sample.C y, SigG.C sig) =>
   y -> T sig t y -> t -> sig t -> sig y -> sig y
multiRelativeZeroPad z ip phase fs x =
   zeroPad multiRelative z ip phase x fs

{-# INLINE multiRelativeConstantPad #-}
multiRelativeConstantPad ::
   (RealField.C t, Sample.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> t -> sig t -> sig y -> sig y
multiRelativeConstantPad ip phase fs x =
   constantPad multiRelative ip phase x fs

{-# INLINE multiRelativeCyclicPad #-}
multiRelativeCyclicPad ::
   (RealField.C t, Sample.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> t -> sig t -> sig y -> sig y
multiRelativeCyclicPad ip phase fs x =
   cyclicPad multiRelative ip phase x fs

{- |
The extrapolation may miss some of the first and some of the last points
-}
{-# INLINE multiRelativeExtrapolationPad #-}
multiRelativeExtrapolationPad ::
   (RealField.C t, Sample.C t, Sample.C y, SigG.C sig) =>
   T sig t y -> t -> sig t -> sig y -> sig y
multiRelativeExtrapolationPad ip phase fs x =
   extrapolationPad multiRelative ip phase x fs
{-
  This example shows pikes, although there shouldn't be any:
   plotList (take 100 $ interpolate (Zero (0::Double)) ipCubic (-0.9::Double) (repeat 0.03) [1,0,1,0.8])
-}

{-* All-in-one interpolation functions -}

{-# INLINE multiRelativeZeroPadConstant #-}
multiRelativeZeroPadConstant ::
   (RealField.C t, Additive.C y, Sample.C t, Sample.C y, SigG.C sig) =>
   t -> sig t -> sig y -> sig y
multiRelativeZeroPadConstant = multiRelativeZeroPad zero constant

{-# INLINE multiRelativeZeroPadLinear #-}
multiRelativeZeroPadLinear ::
   (RealField.C t, Module.C t y, Sample.C t, Sample.C y, SigG.C sig) =>
   t -> sig t -> sig y -> sig y
multiRelativeZeroPadLinear = multiRelativeZeroPad zero linear

{-# INLINE multiRelativeZeroPadCubic #-}
multiRelativeZeroPadCubic ::
   (RealField.C t, Module.C t y, Sample.C t, Sample.C y, SigG.C sig) =>
   t -> sig t -> sig y -> sig y
multiRelativeZeroPadCubic = multiRelativeZeroPad zero cubic


{-* Different kinds of interpolation -}

{-** Hard-wired interpolations -}

data PrefixReader sig a =
   PrefixReader Int (StateT sig Maybe a)

instance Functor (PrefixReader sig) where
   fmap f (PrefixReader count parser) =
      PrefixReader count (fmap f parser)

instance Applicative (PrefixReader sig) where
   pure = PrefixReader 0 . return
   (PrefixReader count0 parser0) <*> (PrefixReader count1 parser1) =
       PrefixReader (count0+count1) (parser0 `ap` parser1)

{-# INLINE getNode #-}
getNode :: (Sample.C y, SigG.C sig) =>
   PrefixReader (sig y) y
getNode = PrefixReader 1 (StateT SigG.viewL)

{-# INLINE fromPrefixReader #-}
fromPrefixReader :: (Sample.C y, SigG.C sig) =>
   String -> Int -> PrefixReader (sig y) (t -> y) -> T sig t y
fromPrefixReader name off (PrefixReader count parser) =
   Cons count off
      (\t xs ->
          maybe
             (error (name ++ " interpolation: not enough nodes"))
             ($t)
             (evalStateT parser xs))

{-| Consider the signal to be piecewise constant. -}
{-# INLINE constant #-}
constant :: (Sample.C y, SigG.C sig) => T sig t y
constant =
   fromPrefixReader "constant" 0 (const <$> getNode)

{-| Consider the signal to be piecewise linear. -}
{-# INLINE linear #-}
linear :: (Module.C t y, Sample.C y, SigG.C sig) => T sig t y
linear =
   fromPrefixReader "linear" 0
      (liftA2
          (\x0 x1 phase -> affineComb phase (x0,x1))
          getNode getNode)

{-| Consider the signal to be piecewise cubic,
    with smooth connections at the nodes.
    It uses a cubic curve which has node values
    x0 at 0 and x1 at 1 and derivatives
    (x1-xm1)/2 and (x2-x0)/2, respectively.
    You can see how it works
    if you evaluate the expression for t=0 and t=1
    as well as the derivative at these points. -}
{-# INLINE cubic #-}
cubic :: (Field.C t, Module.C t y, Sample.C y, SigG.C sig) => T sig t y
cubic =
   fromPrefixReader "cubic" 1
      (liftA4
         (\xm1 x0 x1 x2 t ->
            let lipm12 = affineComb t (xm1,x2)
                lip01  = affineComb t (x0, x1)
                three  = 3 `asTypeOf` t
            in  lip01 + (t*(t-1)/2) *>
                           (lipm12 + (x0+x1) - three *> lip01))
         getNode getNode getNode getNode)

{-# INLINE cubicAlt #-}
cubicAlt :: (Field.C t, Module.C t y, Sample.C y, SigG.C sig) => T sig t y
cubicAlt =
   fromPrefixReader "cubicAlt" 1
      (liftA4
         (\xm1 x0 x1 x2 t ->
          let half = 1/2 `asTypeOf` t
          in  cubicHalf    t  x0 (half *> (x1-xm1)) +
              cubicHalf (1-t) x1 (half *> (x0-x2)))
         getNode getNode getNode getNode)


{- \t -> cubicHalf t x x' has a double zero at 1 and
   at 0 it has value x and steepness x' -}
{-# INLINE cubicHalf #-}
cubicHalf :: (Module.C t y) => t -> y -> y -> y
cubicHalf t x x' =
   (t-1)^2 *> ((1+2*t)*>x + t*>x')



{-** Interpolation based on piecewise defined functions -}

{-# INLINE piecewise #-}
piecewise :: (Module.C t y, Sample.C t, Sample.C y, SigG.C sig) =>
   Int -> [t -> t] -> T sig t y
piecewise center ps =
   Cons (length ps) (center-1)
      (\t -> linearComb (SigG.reverse (SigG.fromList (map ($t) ps))))

{-# INLINE piecewiseConstant #-}
piecewiseConstant ::
   (Module.C t y, Sample.C t, Sample.C y, SigG.C sig) => T sig t y
piecewiseConstant =
   piecewise 1 [const 1]

{-# INLINE piecewiseLinear #-}
piecewiseLinear ::
   (Module.C t y, Sample.C t, Sample.C y, SigG.C sig) => T sig t y
piecewiseLinear =
   piecewise 1 [id, (1-)]

{-# INLINE piecewiseCubic #-}
piecewiseCubic ::
   (Field.C t, Module.C t y, Sample.C t, Sample.C y, SigG.C sig) => T sig t y
piecewiseCubic =
   piecewise 2 $
      Ctrl.cubicFunc (0,(0,0))    (1,(0,1/2)) :
      Ctrl.cubicFunc (0,(0,1/2))  (1,(1,0)) :
      Ctrl.cubicFunc (0,(1,0))    (1,(0,-1/2)) :
      Ctrl.cubicFunc (0,(0,-1/2)) (1,(0,0)) :
      []

{-
GNUPlot.plotList [] $ take 100 $ interpolate (Zero 0) piecewiseCubic (-2.3 :: Double) (repeat 0.1) [2,1,2::Double]
-}


{-** Interpolation based on arbitrary functions -}

{- | with this wrapper you can use the collection of interpolating functions from Donadio's DSP library -}
{-# INLINE function #-}
function :: (Module.C t y, Sample.C t, Sample.C y, SigG.C sig) =>
      (Int,Int)   {- ^ @(left extent, right extent)@, e.g. @(1,1)@ for linear hat -}
   -> (t -> t)
   -> T sig t y
function (left,right) f =
   let len = left+right
   in  Cons len left
          (\t -> linearComb $ SigG.reverse $
               SigG.map
                  (\x -> f (t + fromIntegral x))
                  (SigG.take len (SigG.iterate succ (-left))))
{-
GNUPlot.plotList [] $ take 300 $ interpolate (Zero 0) (function (1,1) (\x -> exp (-6*x*x))) (-2.3 :: Double) (repeat 0.03) [2,1,2::Double]
-}


-- cf. Module.linearComb
{-# INLINE linearComb #-}
linearComb ::
   (SigG.C sig, Sample.C t, Sample.C y, Module.C t y) =>
   sig t -> sig y -> y
linearComb ts ys =
   SigG.sum $ SigG.zipWith (*>) ts ys



{-* Helper functions -}

{-# INLINE delayPad #-}
delayPad :: (Sample.C y, SigG.C sig) => y -> Int -> sig y -> sig y
delayPad z n =
   if n<0 then SigG.drop (negate n) else SigG.append (SigG.replicate n z)