{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Synthesizer.Plain.Control where

import Synthesizer.Plain.Displacement (raise)

import qualified Synthesizer.Plain.Signal as Sig

import qualified Algebra.Module                as Module
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.RealRing              as RealRing
import qualified Algebra.Field                 as Field
import qualified Algebra.Ring                  as Ring
import qualified Algebra.Additive              as Additive

import Algebra.Module((*>))

import Number.Complex (cis,real)
-- import qualified Number.Complex as Complex
import Data.List (zipWith4, tails, )
import Data.List.HT (iterateAssociative, )

import qualified Prelude as P
import NumericPrelude.Base
import NumericPrelude.Numeric


{- * Control curve generation -}

constant :: y -> Sig.T y
constant = repeat


linear :: Additive.C y =>
      y   {-^ steepness -}
   -> y   {-^ initial value -}
   -> Sig.T y {-^ linear progression -}
linear d y0 = iterate (d+) y0

{- |
Minimize rounding errors by reducing number of operations per element
to a logarithmuc number.
-}
linearMultiscale :: Additive.C y =>
      y
   -> y
   -> Sig.T y
linearMultiscale = curveMultiscale (+)

{- |
Linear curve starting at zero.
-}
linearMultiscaleNeutral :: Additive.C y =>
      y
   -> Sig.T y
linearMultiscaleNeutral slope =
   curveMultiscaleNeutral (+) slope zero

{- |
As stable as the addition of time values.
-}
linearStable :: Ring.C y =>
      y
   -> y
   -> Sig.T y
linearStable d y0 =
   curveStable (d*) (+) 1 y0


{- |
It computes the same like 'linear' but in a numerically more stable manner,
namely using a subdivision scheme.
The division needed is a division by two.

0       4       8
0   2   4   6   8
0 1 2 3 4 5 6 7 8
-}
linearMean :: Field.C y =>
      y
   -> y
   -> Sig.T y
linearMean d y0 = y0 :
   foldr (\pow xs -> y0+pow : linearSubdivision xs)
         unreachable (iterate (2*) d)

{- | Intersperse linearly interpolated values. -}
linearSubdivision :: Field.C y =>
      Sig.T y
   -> Sig.T y
linearSubdivision = subdivide (\x0 x1 -> (x0+x1)/2)


{- |
Linear curve of a fixed length.
The final value is not actually reached,
instead we stop one step before.
This way we can concatenate several lines
without duplicate adjacent values.
-}
line :: Field.C y =>
      Int     {-^ length -}
   -> (y,y)   {-^ initial and final value -}
   -> Sig.T y {-^ linear progression -}
line n (y0,y1) =
   take n $ linear ((y1-y0) / fromIntegral n) y0



exponential, exponentialMultiscale, exponentialStable :: Trans.C y =>
      y   {-^ time where the function reaches 1\/e of the initial value -}
   -> y   {-^ initial value -}
   -> Sig.T y {-^ exponential decay -}
exponential time = iterate (* exp (- recip time))
exponentialMultiscale time = curveMultiscale (*) (exp (- recip time))
exponentialStable time = exponentialStableGen exp (- recip time)

exponentialMultiscaleNeutral :: Trans.C y =>
      y   {-^ time where the function reaches 1\/e of the initial value -}
   -> Sig.T y {-^ exponential decay -}
exponentialMultiscaleNeutral time =
   curveMultiscaleNeutral (*) (exp (- recip time)) one

exponential2, exponential2Multiscale, exponential2Stable :: Trans.C y =>
      y   {-^ half life -}
   -> y   {-^ initial value -}
   -> Sig.T y {-^ exponential decay -}
exponential2 halfLife = iterate (*  0.5 ** recip halfLife)
exponential2Multiscale halfLife = curveMultiscale (*) (0.5 ** recip halfLife)
exponential2Stable halfLife = exponentialStableGen (0.5 **) (recip halfLife)

exponential2MultiscaleNeutral :: Trans.C y =>
      y   {-^ half life -}
   -> Sig.T y {-^ exponential decay -}
exponential2MultiscaleNeutral halfLife =
   curveMultiscaleNeutral (*) (0.5 ** recip halfLife) one


exponentialFromTo, exponentialFromToMultiscale :: Trans.C y =>
      y   {-^ time where the function reaches 1\/e of the initial value -}
   -> y   {-^ initial value -}
   -> y   {-^ value after given time -}
   -> Sig.T y {-^ exponential decay -}
exponentialFromTo time y0 y1 =
   iterate (*  (y1/y0) ** recip time) y0
exponentialFromToMultiscale time y0 y1 =
   curveMultiscale (*) ((y1/y0) ** recip time) y0


exponentialStableGen :: (Ring.C y, Ring.C t) =>
      (t -> y)
   -> t
   -> y
   -> Sig.T y
exponentialStableGen expFunc = curveStable expFunc (*)




{-| This is an extension of 'exponential' to vectors
    which is straight-forward but requires more explicit signatures.
    But since it is needed rarely I setup a separate function. -}
vectorExponential :: (Trans.C y, Module.C y v) =>
       y  {-^ time where the function reaches 1\/e of the initial value -}
   ->  v  {-^ initial value -}
   -> Sig.T v {-^ exponential decay -}
vectorExponential time y0 = iterate (exp (-1/time) *>) y0

vectorExponential2 :: (Trans.C y, Module.C y v) =>
       y  {-^ half life -}
   ->  v  {-^ initial value -}
   -> Sig.T v {-^ exponential decay -}
vectorExponential2 halfLife y0 = iterate (0.5**(1/halfLife) *>) y0



cosine, cosineMultiscale, cosineSubdiv, cosineStable :: Trans.C y =>
       y  {-^ time t0 where  1 is approached -}
   ->  y  {-^ time t1 where -1 is approached -}
   -> Sig.T y {-^ a cosine wave where one half wave is between t0 and t1 -}
cosine = cosineWithSlope $
   \d x -> map cos (linear d x)

cosineMultiscale = cosineWithSlope $
   \d x -> map real (curveMultiscale (*) (cis d) (cis x))


{-
  cos (a-b) = cos a * cos b + sin a * sin b
  cos (a+b) = cos a * cos b - sin a * sin b
  cos  a    = (cos (a-b) + cos (a+b)) / (2 * cos b)

  Problem: (cos b) might be close to zero,
  example: Syn.cosineStable 1 (9::Double)
-}
cosineSubdiv =
   let aux d y0 =
          cos y0 :
            foldr (\pow xs -> cos(y0+pow) : cosineSubdivision pow xs)
                  unreachable (iterate (2*) d)
   in  cosineWithSlope aux

cosineSubdivision :: Trans.C y =>
      y
   -> Sig.T y
   -> Sig.T y
cosineSubdivision angle =
   let k = recip (2 * cos angle)
   in  subdivide (\x0 x1 -> (x0+x1)*k)

cosineStable = cosineWithSlope $
   \d x -> map real (exponentialStableGen cis d (cis x))


cosineWithSlope :: Trans.C y =>
      (y -> y -> signal)
   ->  y
   ->  y
   -> signal
cosineWithSlope c t0 t1 =
   let inc = pi/(t1-t0)
   in  c inc (-t0*inc)


cubicHermite :: Field.C y => (y, (y,y)) -> (y, (y,y)) -> Sig.T y
cubicHermite node0 node1 =
   map (cubicFunc node0 node1) (linear 1 0)

{- |
0                                     16
0               8                     16
0       4       8         12          16
0   2   4   6   8   10    12    14    16
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
-}
cubicFunc :: Field.C y => (y, (y,y)) -> (y, (y,y)) -> y -> y
cubicFunc (t0, (y0,dy0)) (t1, (y1,dy1)) t =
   let dt  = t0-t1
       dt0 = t-t0
       dt1 = t-t1
       x0  = dt1^2
       x1  = dt0^2
   in  ((dy0*dt0 + y0 * (1-2/dt*dt0)) * x0 +
        (dy1*dt1 + y1 * (1+2/dt*dt1)) * x1) / dt^2
{-
cubic t0 (y0,dy0) t1 (y1,dy1) t =
   let x0 = ((t-t1) / (t0-t1))^2
       x1 = ((t-t0) / (t1-t0))^2
   in  y0 * x0 + y1 * x1 +
       (dy0 - y0*2/(t0-t1)) * (t-t0)*x0 +
       (dy1 - y1*2/(t1-t0)) * (t-t1)*x1
-}

cubicHermiteStable :: Field.C y => (y, (y,y)) -> (y, (y,y)) -> Sig.T y
cubicHermiteStable node0 node1 =
   cubicFunc node0 node1 0 :
      foldr (\pow xs ->
                cubicFunc node0 node1 pow : head xs :
                cubicFunc node0 node1 (3*pow) : cubicSubdivision xs)
            unreachable (iterate (2*) 1)

cubicSubdivision :: Field.C y => Sig.T y -> Sig.T y
cubicSubdivision xs =
   let xs0:xs1:xs2:xs3:_ = tails xs
       inter = zipWith4 (\x0 x1 x2 x3 -> (9*(x1+x2) - (x0+x3))/16)
                        xs0 xs1 xs2 xs3
   in  head xs1 : flattenPairs (zip inter xs2)

{-
            foldr (\(pow0:pow1:_) ~(_:xs) ->
                      cos (y0+pow0) : cos (y0+pow1) : cos (y0+pow0+pow1) :
                         cosineSubdivision pow0 xs)
                  unreachable (tails (iterate (2*) d))
-}


{-
maybe cubicHermite could also be implemented in a Multiscale manner
using a difference scheme.

cubicHermiteMultiscale :: Field.C y => (y, (y,y)) -> (y, (y,y)) -> Sig.T y
cubicHermiteMultiscale node0@(t0,y0) node1@(t1,y1) =
   let -- could be inlined and simplified
       ys = map (cubicFunc node0 node1) [0,1,2,3]
       (d0:d1:d2:d3:_) = iterate (mapAdjacent substract) ys

I thought multiplying difference schemes could help somehow,
but it doesn't. :-(

cubicHermiteMultiscale

Leibniz rule for differences

D3(s+r) = D0(s)*D3(r) + 3*D1(s)*D2(r) + 3*D2(s)*D1(r) + D3(s)*D0(r)


mulDiffs4 :: Ring.C a => (a,a,a,a) -> (a,a,a,a) -> (a,a,a,a)
mulDiffs4 (r0,r1,r2,r3) (s0,s1,s2,s3) =
   (r0*s0,
    r0*s1 +   r1*s0,
    r0*s2 + 2*r1*s1 +   r2*s0,
    r0*s3 + 3*r1*s2 + 3*r2*s1 + r3*s0)

mulDiffs4zero :: Ring.C a => (a,a,a) -> (a,a,a) -> (a,a,a)
mulDiffs4zero (r0,r1,r2,r3) (s0,s1,s2,s3) =
   (r0*s0,
    r0*s1 +   r1*s0,
    r0*s2 + 2*r1*s1 +   r2*s0,
    r0*s3 + 3*r1*s2 + 3*r2*s1 + r3*s0)

mulDiffs3 :: Ring.C a => (a,a,a) -> (a,a,a) -> (a,a,a)
mulDiffs3 (r0,r1,r2) (s0,s1,s2) =
   (r0*s0,
    r0*s1 +   r1*s0,
    r0*s2 + 2*r1*s1 +   r2*s0)

mulDiffs3Karatsuba :: Ring.C a => (a,a,a) -> (a,a,a) -> (a,a,a)
mulDiffs3Karatsuba (r0,r1,r2) (s0,s1,s2) =
   let r0s0 = r0*s0
       r1s1 = r1*s1
   in  (r0s0,
        (r0+r1)*(s0+s1) - r0s0 - r1s1,
        r0*s2 + 2*r1s1 + r2*s0)
-}



{-# DEPRECATED Control "use Synthesizer.State.Piece instead" #-}
{- |
The curve type of a piece of a piecewise defined control curve.
-}
data Control y =
     CtrlStep
   | CtrlLin
   | CtrlExp {ctrlExpSaturation :: y}
   | CtrlCos
   | CtrlCubic {ctrlCubicGradient0 :: y,
                ctrlCubicGradient1 :: y}
   deriving (Eq, Show)

{- |
The full description of a control curve piece.
-}
data ControlPiece y =
     ControlPiece {pieceType :: Control y,
                   pieceY0 :: y,
                   pieceY1 :: y,
                   pieceDur :: y}
   deriving (Eq, Show)


newtype PieceRightSingle y = PRS y
newtype PieceRightDouble y = PRD y

type ControlDist y = (y, Control y, y)


-- precedence and associativity like (:)
infixr 5 -|#, #|-, =|#, #|=, |#, #|

{- |
The 6 operators simplify constructing a list of @ControlPiece a@.
The description consists of nodes (namely the curve values at nodes)
and the connecting curve types.
The naming scheme is as follows:
In the middle there is a bar @|@.
With respect to the bar,
the pad symbol @\#@ is at the side of the curve type,
at the other side there is nothing, a minus sign @-@, or an equality sign @=@.

 (1) Nothing means that here is the start or the end node of a curve.

 (2) Minus means that here is a node where left and right curve meet at the same value.
     The node description is thus one value.

 (3) Equality sign means that here is a split node,
     where left and right curve might have different ending and beginning values, respectively.
     The node description consists of a pair of values.
-}

-- the leading space is necessary for the Haddock parser

( #|-) :: (y, Control y) -> (PieceRightSingle y, [ControlPiece y]) ->
   (ControlDist y, [ControlPiece y])
(d,c) #|- (PRS y1, xs)  =  ((d,c,y1), xs)

(-|#) :: y -> (ControlDist y, [ControlPiece y]) ->
   (PieceRightSingle y, [ControlPiece y])
y0 -|# ((d,c,y1), xs)  =  (PRS y0, ControlPiece c y0 y1 d : xs)

( #|=) :: (y, Control y) -> (PieceRightDouble y, [ControlPiece y]) ->
   (ControlDist y, [ControlPiece y])
(d,c) #|= (PRD y1, xs)  =  ((d,c,y1), xs)

(=|#) :: (y,y) -> (ControlDist y, [ControlPiece y]) ->
   (PieceRightDouble y, [ControlPiece y])
(y01,y10) =|# ((d,c,y11), xs)  =  (PRD y01, ControlPiece c y10 y11 d : xs)

( #|) :: (y, Control y) -> y ->
   (ControlDist y, [ControlPiece y])
(d,c) #| y1  =  ((d,c,y1), [])

(|#) :: y -> (ControlDist y, [ControlPiece y]) ->
   [ControlPiece y]
y0 |# ((d,c,y1), xs)  =  ControlPiece c y0 y1 d : xs


piecewise :: (Trans.C y, RealRing.C y) =>
   [ControlPiece y] -> Sig.T y
piecewise xs =
   let ts = scanl (\(_,fr) d -> splitFraction (fr+d))
                  (0,1) (map pieceDur xs)
   in  concat (zipWith3
          (\n t (ControlPiece c yi0 yi1 d) ->
               piecewisePart yi0 yi1 t d n c)
          (map fst (tail ts)) (map (subtract 1 . snd) ts)
          xs)


piecewisePart :: (Trans.C y) =>
   y -> y -> y -> y -> Int -> Control y -> Sig.T y
piecewisePart y0 y1 t0 d n ctrl =
   take n
      (case ctrl of
         CtrlStep  -> constant y0
         CtrlLin   -> let s = (y1-y0)/d in linearStable s (y0-t0*s)
         CtrlExp s -> let y0' = y0-s; y1' = y1-s; yd = y0'/y1'
                      in  raise s (exponentialStable (d / log yd)
                                           (y0' * yd**(t0/d)))
         CtrlCos   -> map (\y -> (1+y)*(y0/2)+(1-y)*(y1/2))
                          (cosineStable t0 (t0+d))
         CtrlCubic yd0 yd1 ->
            cubicHermiteStable (t0,(y0,yd0)) (t0+d,(y1,yd1)))

{-
  exp (-1/time) == yd**(-1/d)
  1/time == log yd / d
  time   == d / log yd
-}

{-
  piecewise (0 |# (10.21, CtrlExp 1.1) #|- 1 -|# (10,CtrlExp 0.49) #|- 0.5 -|# (30, CtrlLin) #|- 0.5 -|# (20, CtrlCos) #| 0)

  piecewise (0 |# (10.21, CtrlExp 1.1) #|- 1 -|# (10,CtrlCubic (-0.1) 0) #|- 0.5 -|# (30, CtrlLin) #|- 0.5 -|# (20, CtrlCos) #| 0)
-}


{- * Auxiliary functions -}

curveStable :: (Additive.C t) =>
      (t -> y)
   -> (y -> y -> y)
   -> t
   -> y
   -> Sig.T y
curveStable expFunc op time y0 =
   y0 : map (op y0)
      (foldr
         (\e xs ->
            let k = expFunc e
            in  k : concatMapPair (\x -> (x, op x k)) xs)
       unreachable (iterate double time))

unreachable :: a
unreachable = error "only reachable in infinity"

double :: Additive.C t => t -> t
double t = t+t

concatMapPair :: (a -> (b,b)) -> Sig.T a -> Sig.T b
concatMapPair f = flattenPairs . map f

flattenPairs :: Sig.T (a,a) -> Sig.T a
flattenPairs = foldr (\(a,b) xs -> a:b:xs) []

subdivide :: (y -> y -> y) -> Sig.T y -> Sig.T y
subdivide f xs0@(x:xs1) =
   x : flattenPairs (zipWith (\x0 x1 -> (f x0 x1, x1)) xs0 xs1)
subdivide _ [] = []


concatMapPair' :: (a -> (b,b)) -> Sig.T a -> Sig.T b
concatMapPair' f = concatMap ((\(x,y) -> [x,y]) . f)


curveMultiscale :: (y -> y -> y) -> y -> y -> Sig.T y
curveMultiscale op d y0 =
   y0 : map (op y0) (iterateAssociative op d)


curveMultiscaleNeutral :: (y -> y -> y) -> y -> y -> Sig.T y
curveMultiscaleNeutral op d neutral =
   neutral : iterateAssociative op d