{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-}
module Synthesizer.Generic.Control where

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

import Synthesizer.Generic.Displacement (raise)

import qualified Algebra.Module                as Module
import qualified Algebra.Transcendental        as Trans
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.Module((*>))

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

import qualified Prelude as P
import PreludeBase
import NumericPrelude


{- * Control curve generation -}

constant :: (Sample.C y, SigG.C sig) => y -> sig y
constant = SigG.repeat


linear :: (Additive.C y, Sample.C y, SigG.C sig) =>
      y   {-^ steepness -}
   -> y   {-^ initial value -}
   -> sig y {-^ linear progression -}
linear d y0 = SigG.iterate (d+) y0

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

{- |
Linear curve starting at zero.
-}
linearMultiscaleNeutral :: (Additive.C y, Sample.C y, SigG.C sig) =>
      y
   -> sig y
linearMultiscaleNeutral slope =
   curveMultiscaleNeutral (+) slope zero


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

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

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

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




{-| 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, Sample.C v, SigG.C sig) =>
       y  {-^ time where the function reaches 1\/e of the initial value -}
   ->  v  {-^ initial value -}
   -> sig v {-^ exponential decay -}
vectorExponential time y0 = SigG.iterate (exp (-1/time) *>) y0

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



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

cosineMultiscale = cosineWithSlope $
   \d x -> SigG.map real (curveMultiscale (*) (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, Sample.C y, SigG.C sig) =>
   (y, (y,y)) -> (y, (y,y)) -> sig y
cubicHermite node0 node1 =
   SigG.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
-}



{- |
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, RealField.C y, Sample.C y, SigG.C sig) =>
   [ControlPiece y] -> sig y
piecewise xs =
   let ts = scanl (\(_,fr) d -> splitFraction (fr+d))
                  (0,1) (map pieceDur xs)
   in  SigG.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, Sample.C y, SigG.C sig) =>
   y -> y -> y -> y -> Int -> Control y -> sig y
piecewisePart y0 y1 t0 d n ctrl =
   SigG.take n
      (case ctrl of
         CtrlStep  -> constant y0
         CtrlLin   -> let s = (y1-y0)/d in linearMultiscale s (y0-t0*s)
         CtrlExp s -> let y0' = y0-s; y1' = y1-s; yd = y0'/y1'
                      in  raise s (exponentialMultiscale (d / log yd)
                                           (y0' * yd**(t0/d)))
         CtrlCos   -> SigG.map
                          (\y -> (1+y)*(y0/2)+(1-y)*(y1/2))
                          (cosineMultiscale t0 (t0+d))
         CtrlCubic yd0 yd1 ->
            cubicHermite (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 -}


curveMultiscale :: (Sample.C y, SigG.C sig) =>
   (y -> y -> y) -> y -> y -> sig y
curveMultiscale op d y0 =
   SigG.cons y0 (SigG.map (op y0) (SigG.iterateAssoc op d))


curveMultiscaleNeutral :: (Sample.C y, SigG.C sig) =>
   (y -> y -> y) -> y -> y -> sig y
curveMultiscaleNeutral op d neutral =
   SigG.cons neutral (SigG.iterateAssoc op d)