```{-# 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.Plain.Interpolation where

import qualified Synthesizer.Plain.Control as Ctrl
import qualified Synthesizer.Plain.Signal  as Sig
import qualified Synthesizer.Plain.Filter.NonRecursive as FiltNR

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 Algebra.Module((*>))
import Data.Maybe (fromMaybe)
import Synthesizer.Utility (viewListL, viewListR, affineComb, )
import Synthesizer.ApplicativeUtility (liftA4, )

import Control.Monad.State (StateT(StateT), evalStateT, replicateM_, ap, guard, )
import Control.Applicative (Applicative(pure, (<*>)), (<\$>), liftA2, )

import PreludeBase
import NumericPrelude

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

{-* Interpolation with various padding methods -}

(T t y -> t -> Sig.T y -> a) ->
y -> T t y -> t -> Sig.T y -> a
zeroPad interpolate z ip phase x =
let (phInt, phFrac) = splitFraction phase
in  interpolate ip phFrac
(FiltNR.delayPad z (offset ip - phInt) (x ++ repeat z))

(T t y -> t -> Sig.T y -> a) ->
T t y -> t -> Sig.T y -> a
constantPad interpolate ip phase x =
let (phInt, phFrac) = splitFraction phase
do (xFirst,_) <- viewListL x
(xBody,xLast) <- viewListR x
return (FiltNR.delayPad xFirst (offset ip - phInt) (xBody ++ repeat xLast))
in  interpolate ip phFrac

{- |
Only for finite input signals.
-}
(T t y -> t -> Sig.T y -> a) ->
T t y -> t -> Sig.T y -> a
cyclicPad interpolate ip phase x =
let (phInt, phFrac) = splitFraction phase
in  interpolate ip phFrac
(drop (mod (phInt - offset ip) (length x)) (cycle x))

{- |
The extrapolation may miss some of the first and some of the last points
-}
(T t y -> t -> Sig.T y -> a) ->
T t y -> t -> Sig.T y -> a
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 -}

skip :: (RealField.C t) =>
T t y -> (t, Sig.T y) -> (t, Sig.T y)
skip ip (phase0, x0) =
let (n, frac) = splitFraction phase0
(m, x1) = Sig.dropMarginRem (number ip) n x0
in  (fromIntegral m + frac, x1)

single :: (RealField.C t) =>
T t y -> t -> Sig.T 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])
-}

-- | alternative implementation of 'single'
singleRec :: (Ord t, Ring.C t) =>
T t y -> t -> Sig.T y -> y
singleRec ip phase x =
-- check if we are leaving the current interval
maybe
(func ip phase x)
(singleRec ip (phase - 1))
(do (_,xs) <- viewListL x
guard (phase >= 1 && minLength (number ip) xs)
return xs)

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

{- | All values of frequency control must be non-negative. -}
multiRelative :: (RealField.C t) =>
T t y -> t -> Sig.T y -> Sig.T t -> Sig.T y
multiRelative ip phase0 x0 =
map (uncurry (func ip)) .
scanl
(\(phase,x) freq -> skip ip (phase + freq, x))
(skip ip (phase0,x0))

y -> T t y -> t -> Sig.T t -> Sig.T y -> Sig.T y
multiRelativeZeroPad z ip phase fs x =
zeroPad multiRelative z ip phase x fs

T t y -> t -> Sig.T t -> Sig.T y -> Sig.T y
multiRelativeConstantPad ip phase fs x =
constantPad multiRelative ip phase x fs

T t y -> t -> Sig.T t -> Sig.T y -> Sig.T 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
-}
T t y -> t -> Sig.T t -> Sig.T y -> Sig.T 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 -}

(RealField.C t, Additive.C y) => t -> Sig.T t -> Sig.T y -> Sig.T y

(RealField.C t, Module.C t y) => t -> Sig.T t -> Sig.T y -> Sig.T y

(RealField.C t, Module.C t y) => t -> Sig.T t -> Sig.T y -> Sig.T y

{-* Different kinds of interpolation -}

{-** Hard-wired interpolations -}

PrefixReader Int (StateT (Sig.T y) Maybe a)

fmap f (PrefixReader count parser) =

-- this is a MonadWriter with Sum monoid
pure = PrefixReader 0 . return

getNode = PrefixReader 1 (StateT viewListL)

fromPrefixReader :: String -> Int -> PrefixReader y (t -> y) -> T t y
Cons count off
(\t xs ->
maybe
(error (name ++ " interpolation: not enough nodes"))
(\$t)
(evalStateT parser xs))

{-| Consider the signal to be piecewise constant. -}
constant :: T t y
constant =
fromPrefixReader "constant" 0 (const <\$> getNode)

{-| Consider the signal to be piecewise linear. -}
linear :: (Module.C t y) => T t y
linear =
(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. -}
cubic :: (Field.C t, Module.C t y) => T t y
cubic =
(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)

cubicAlt :: (Field.C t, Module.C t y) => T t y
cubicAlt =
(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' -}
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 -}

piecewise :: (Module.C t y) =>
Int -> [t -> t] -> T t y
piecewise center ps =
Cons (length ps) (center-1)
(\t -> Module.linearComb (reverse (map (\$t) ps)))

piecewiseConstant :: (Module.C t y) => T t y
piecewiseConstant =
piecewise 1 [const 1]

piecewiseLinear :: (Module.C t y) => T t y
piecewiseLinear =
piecewise 1 [id, (1-)]

piecewiseCubic :: (Field.C t, Module.C t y) => T 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 -}
function :: (Module.C t y) =>
(Int,Int)   {- ^ @(left extent, right extent)@, e.g. @(1,1)@ for linear hat -}
-> (t -> t)
-> T t y
function (left,right) f =
let len = left+right
in  Cons len left
(\t -> Module.linearComb (reverse
(map (\x -> f (t + fromIntegral x)) (take len [(-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]
-}

{-* Helper functions -}

{-| Test if a list has at least @n@ elements
make sure that @n@ is non-negative -}
minLength :: Int -> Sig.T y -> Bool
minLength n xs =
maybe False (const True) (evalStateT (replicateM_ n (StateT viewListL)) xs)
```