```{-# LANGUAGE NoImplicitPrelude #-}
{- |
Special interpolations defined in terms of Module operations.
-}
module Synthesizer.Interpolation.Module (
T,
constant,
linear,
cubic,
cubicAlt,
piecewise,
piecewiseConstant,
piecewiseLinear,
piecewiseCubic,
function,
) where

import qualified Synthesizer.State.Signal  as Sig
import qualified Synthesizer.Plain.Control as Ctrl

import Synthesizer.Interpolation (
constant,
)

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 Control.Applicative (liftA2, )
import Synthesizer.ApplicativeUtility (liftA4, )
import Synthesizer.Utility (affineComb, )

import PreludeBase
import NumericPrelude

{-| Consider the signal to be piecewise linear. -}
{-# INLINE 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.
-}
{-# INLINE cubic #-}
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)

{- |
The interpolators for module operations
do not simply compute a straight linear combination of some vectors.
This is efficient whenever scaling and addition is cheap.
In this case they might save multiplications.
I can't say much about numeric cancellations, however.
-}
{-# INLINE cubicAlt #-}
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 slope 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) =>
Int -> [t -> t] -> T t y
piecewise center ps =
cons (length ps) (center-1)
(\t -> Sig.linearComb (Sig.fromList (map (\$t) (reverse ps))))

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

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

{-# INLINE piecewiseCubic #-}
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 -}
{-# INLINE function #-}
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
ps  = Sig.take len \$ Sig.iterate pred (pred right)
-- ps = Sig.reverse \$ Sig.take len \$ Sig.iterate succ (-left)
in  cons len left
(\t -> Sig.linearComb \$
Sig.map (\x -> f (t + fromIntegral x)) ps)
{-
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]
-}
```