{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Synthesizer.Plain.Control (
   constant,
   linear,
   linearMultiscale,
   linearMultiscaleNeutral,
   linearStable,
   linearMean,
   line,
   exponential, exponentialMultiscale, exponentialStable,
   exponentialMultiscaleNeutral,
   exponential2, exponential2Multiscale, exponential2Stable,
   exponential2MultiscaleNeutral,
   exponentialFromTo, exponentialFromToMultiscale,
   vectorExponential,
   vectorExponential2,
   cosine, cosineMultiscale, cosineSubdiv, cosineStable,
   cubicHermite,
   cubicHermiteStable,

   -- used in Analysis
   curveMultiscale,
   curveMultiscaleNeutral,
   -- used in Generic.Control, Interpolation.Module
   cubicFunc,
   cosineWithSlope,
   ) where

import qualified Synthesizer.Plain.Signal as Sig

import Data.List (zipWith4, tails, )
import Data.List.HT (iterateAssociative, )

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

import Number.Complex (cis,real, )

import NumericPrelude.Numeric
import NumericPrelude.Base


{- * 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)
-}



{- * 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