{-# OPTIONS -O2 -fno-implicit-prelude #-}
{- |
Copyright   :  (c) Henning Thielemann 2006
License     :  GPL

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes


Avoid importing this module.
Better use functions from
"Synthesizer.Plain.Oscillator" and
"Synthesizer.Basic.Wave"

Input data is interpreted as samples of data on a cylinder
in the following form:

> |*          |
> |   *       |
> |      *    |
> |         * |
> | *         |
> |    *      |
> |       *   |
> |          *|
> |  *        |
> |     *     |
> |        *  |


> -----------
> *
>     *
>         *
>  *
>      *
>          *
>   *
>       *
>           *
>    *
>        *
> -----------

We have to interpolate in the parallelograms.

-}
module Synthesizer.Plain.ToneModulation where

import qualified Synthesizer.Basic.Phase as Phase

import qualified Synthesizer.Plain.Interpolation as Interpolation
-- import qualified Data.Array as Array
import Data.Array (Array, (!), listArray, )

-- import qualified Algebra.Transcendental        as Trans
import qualified Algebra.RealField             as RealField
import qualified Algebra.Field                 as Field
-- import qualified Algebra.Real                  as Real
import qualified Algebra.Ring                  as Ring
import qualified Algebra.Additive              as Additive

import qualified Number.NonNegative       as NonNeg
import qualified Number.NonNegativeChunky as Chunky

import Synthesizer.Utility (viewListL, viewListR, clip, mapPair, )
import Control.Monad (guard)

import qualified Data.List as List

import NumericPrelude.List (replicateMatch, takeMatch, )
import NumericPrelude

-- import qualified Prelude as P
import PreludeBase



-- * general helpers

interpolateCell ::
   Interpolation.T a y ->
   Interpolation.T b y ->
   (a, b) ->
   [[y]] -> y
interpolateCell ipLeap ipStep (qLeap,qStep) =
   Interpolation.func ipStep qStep .
   map (Interpolation.func ipLeap qLeap)


{- |
Convert from the (shape,phase) parameter pair
to the index within a wave (step) and the index of a wave (leap)
in the sampled prototype tone.
-}
untangleShapePhase :: (Field.C a) =>
   Int -> a -> (a, a) -> (a, a)
untangleShapePhase periodInt period (shape,phase) =
   let leap = shape/period - phase
       step = shape - leap * fromIntegral periodInt
   in  (leap, step)

untangleShapePhaseAnalytic :: (Field.C a) =>
   Int -> a -> (a, a) -> (a, a)
untangleShapePhaseAnalytic periodInt period (shape,phase) =
   let periodRound = fromIntegral periodInt
       vLeap = (periodRound, periodRound-period)
       vStep = (1,1)
   in  solveSLE2 (vLeap,vStep) (shape,period*phase)

{-
Cramer's rule

see HTam/Numerics/ZeroFinder/Root, however the matrix is transposed
-}
solveSLE2 :: Field.C a => ((a,a), (a,a)) -> (a,a) -> (a,a)
solveSLE2 a@(a0,a1) b =
   let det = det2 a
   in  (det2 (b, a1) / det,
        det2 (a0, b) / det)

det2 :: Ring.C a => ((a,a), (a,a)) -> a
det2 ((a00,a10),(a01,a11)) =
   a00*a11 - a10*a01

{-
transpose :: ((a,a), (a,a)) -> ((a,a), (a,a))
transpose ((a00,a10),(a01,a11)) = ((a00,a01),(a10,a11))
-}


flattenShapePhase :: RealField.C a =>
      Int
   -> a
   -> (a, Phase.T a)
   -> (Int, (a, a))
flattenShapePhase periodInt period (shape,phase) =
   let (xShape,xWave) =
          untangleShapePhase periodInt period (shape, Phase.toRepresentative phase)
       (nLeap,qLeap) = splitFraction xShape
       (nStep,qStep) = splitFraction xWave
       {- reverse solveSLE2 for the shape parameter
          with respect to the rounded (wave,shape) coordinates -}
       n = nStep + nLeap * periodInt
   in  (n,(qLeap,qStep))

shapeLimits :: Ring.C t =>
      Interpolation.T a v
   -> Interpolation.T a v
   -> Int
   -> t
   -> (t, t)
shapeLimits ipLeap ipStep periodInt len =
   let minShape =
          fromIntegral $
          interpolationOffset ipLeap ipStep periodInt +
          periodInt
       maxShape =
          minShape + len -
             fromIntegral
                (Interpolation.number ipStep +
                 Interpolation.number ipLeap * periodInt)
   in  (minShape, maxShape)


interpolationOffset ::
      Interpolation.T a v
   -> Interpolation.T a v
   -> Int
   -> Int
interpolationOffset ipLeap ipStep periodInt =
   Interpolation.offset ipStep +
   Interpolation.offset ipLeap * periodInt




-- * array based shape variable wave

data Prototype a v =
   Prototype {
      protoIpLeap,
      protoIpStep      :: Interpolation.T a v,
      protoIpOffset    :: Int,
      protoPeriod      :: a,
      protoPeriodInt   :: Int,
      protoShapeLimits :: (a,a),
      protoArray       :: Array Int v
   }


makePrototype :: (RealField.C a) =>
   Interpolation.T a v ->
   Interpolation.T a v ->
   a -> [v] -> Prototype a v
makePrototype ipLeap ipStep period tone =
   let periodInt = round period
       ipOffset =
          interpolationOffset ipLeap ipStep periodInt
       len = length tone
       (lower,upper) =
          shapeLimits ipLeap ipStep periodInt len
       limits =
          if lower > upper
            then error "min>max"
            else
              (fromIntegral lower, fromIntegral upper)

       arr = listArray (0, pred len) tone

   in  Prototype {
          protoIpLeap      = ipLeap,
          protoIpStep      = ipStep,
          protoIpOffset    = ipOffset,
          protoPeriod      = period,
          protoPeriodInt   = periodInt,
          protoShapeLimits = limits,
          protoArray       = arr
       }

sampledToneCell :: (RealField.C a) =>
   Prototype a v -> a -> Phase.T a -> ((a,a),[[v]])
sampledToneCell p shape phase =
   let (n, q) =
          flattenShapePhase (protoPeriodInt p) (protoPeriod p)
             (uncurry clip (protoShapeLimits p) shape, phase)
   in  (q,
        map (map (protoArray p ! ) . iterate (protoPeriodInt p +)) $
        enumFrom (n - protoIpOffset p))

sampledToneAltCell :: (RealField.C a) =>
   Prototype a v -> a -> Phase.T a -> ((a,a),[[v]])
sampledToneAltCell p shape phase =
   let (n, q) =
          flattenShapePhase (protoPeriodInt p) (protoPeriod p)
             (uncurry clip (protoShapeLimits p) shape, phase)
   in  (q,
        iterate (drop (protoPeriodInt p)) $
        map (protoArray p ! ) (enumFrom (n - protoIpOffset p)))


{-
  M = ((1,1)^T, (periodRound, period-periodRound)^T)

  equation for the line
   0 = (nStep  - offset ipStep) +
       (nLeap - offset ipLeap) * periodInt

   <(1,periodInt), (offset ipStep, offset ipLeap)>
        = <(1,periodInt), (nStep,nLeap)>
   d = <a,x>
     = <a,M^-1*M*x>
     = <(M^-T)*a,M*x>
     = <(M^-T)*a,y>
   b = (M^-T)*a
   required:
      y0 such that y1=0
      y0 such that y1=period

   The line {x : d = <a,x>} converted to (shape,phase) coordinates
   has constant shape and meets all phases.
-}



-- * lazy oscillator


oscillatorCells :: (RealField.C t) =>
    Interpolation.T t y ->
    Interpolation.T t y ->
    t -> [y] -> (t,[t]) -> (Phase.T t,[t]) -> [((t,t),[[y]])]
oscillatorCells
       ipLeap ipStep period sampledTone shapes freqs =
    let periodInt = round period
        ptrs =
           List.transpose $
           takeWhile (not . null) $
           iterate (drop periodInt) sampledTone
        ipOffset =
           interpolationOffset ipLeap ipStep periodInt
        (skip:skips,coords) =
           unzip $
           oscillatorCoords periodInt period
              (limitRelativeShapes ipLeap ipStep periodInt sampledTone shapes)
              freqs
    in  zipWith
           -- n will be zero within the data, it's only needed for extrapolation
           (\(k,q) (n,ptr) ->
             if n>0
               then error "ToneModulation.oscillatorCells: limit of shape parameter is buggy"
               else
                 (q, drop (periodInt+k) ptr))
           coords $
        tail $
        scanl
           {- since we clip the coordinates before calling oscillatorCells
              we do not need 'dropRem', since 'drop' would never go beyond the list end -}
           (\ (n,ptr0) d0 -> dropRem (n+d0) ptr0)
           (0,ptrs)
           ((skip - ipOffset - periodInt) : skips)

dropFrac :: RealField.C i => i -> [a] -> (Int, i, [a])
dropFrac =
   let recurse acc n xt =
          if n>=1
            then
               case xt of
                  _:xs -> recurse (succ acc) (n-1) xs
                  [] -> (acc, n, [])
            else (acc,n,xt)
   in  recurse 0

dropFrac' :: RealField.C i => i -> [a] -> (Int, i, [a])
dropFrac' =
   let recurse acc n xt =
          maybe
             (acc,n,xt)
             (recurse (succ acc) (n-1) . snd)
             (guard (n>=1) >> viewListL xt)
   in  recurse 0

propDropFrac :: (RealField.C i, Eq a) => i -> [a] -> Bool
propDropFrac n xs =
   dropFrac n xs == dropFrac' n xs



dropRem :: Int -> [a] -> (Int, [a])
dropRem =
   let recurse n xt =
          if n>0
            then
               case xt of
                  _:xs -> recurse (pred n) xs
                  [] -> (n, [])
            else (n,xt)
   in  recurse

dropRem' :: Int -> [a] -> (Int, [a])
dropRem' =
   let recurse n xt =
          maybe
             (n,xt)
             (recurse (pred n) . snd)
             (guard (n>0) >> viewListL xt)
   in  recurse

propDropRem :: (Eq a) => Int -> [a] -> Bool
propDropRem n xs =
   dropRem n xs == dropRem' n xs

{-
*Synthesizer.Plain.ToneModulation> Test.QuickCheck.quickCheck (\n xs -> propDropRem n (xs::[Int]))
OK, passed 100 tests.
*Synthesizer.Plain.ToneModulation> Test.QuickCheck.quickCheck (\n xs -> propDropFrac (n::Rational) (xs::[Int]))
OK, passed 100 tests.
-}


oscillatorCoords :: (RealField.C t) =>
    Int -> t -> (t,[t]) -> (Phase.T t, [t]) -> [(Int,(Int,(t,t)))]
oscillatorCoords periodInt period
       (shape0, shapes) (phase, freqs) =
    let shapeOffsets =
           scanl
              (\(_,s) c -> splitFraction (s+c))
              (splitFraction shape0) shapes
        phases =
           let (s:ss) = map (\(n,_) -> fromIntegral n / period) shapeOffsets
           in  freqToPhase
                  (Phase.increment (-s) phase)  -- phase - s
                  (zipWith (-) freqs ss)
    in  zipWith
--           (\(d,s) p -> (d, (s,p)))
           (\(d,s) p -> (d, flattenShapePhase periodInt period (s,p)))
           shapeOffsets
           phases
{-
mapM print $ take 30 $ let period = 1/0.07::Double in oscillatorCoords (round period) period 0 0 (repeat 0.1) (repeat 0.01)

*Synthesizer.Plain.Oscillator> mapM print $ take 30 $ let period = 1/0.07::Rational in oscillatorCoords (round period) period 0 0 (repeat 1) (repeat 0.07)

*Synthesizer.Plain.Oscillator> mapM print $ take 30 $ let period = 1/0.07::Rational in oscillatorCoords (round period) period 0 0 (repeat 0.25) (repeat 0.0175)
-}


-- this function fits better in the Oscillator module
{- |
Convert a list of phase steps into a list of momentum phases
phase is a number in the interval [0,1)
freq contains the phase steps
-}
freqToPhase :: RealField.C a => Phase.T a -> [a] -> [Phase.T a]
freqToPhase phase freq = scanl (flip Phase.increment) phase freq



limitRelativeShapes :: (RealField.C t) =>
    Interpolation.T t y ->
    Interpolation.T t y ->
    Int -> [y] -> (t,[t]) -> (t,[t])
limitRelativeShapes ipLeap ipStep periodInt sampledTone =
    let -- len = List.genericLength sampledTone
        len = Chunky.fromChunks (replicateMatch sampledTone one)
        (minShape, maxShape) = shapeLimits ipLeap ipStep periodInt len
        fromChunky = NonNeg.toNumber   . Chunky.toNumber
        toChunky   = Chunky.fromNumber . NonNeg.fromNumber
    in  mapPair (fromChunky, map fromChunky) .
        uncurry (limitMaxRelativeValuesNonNeg maxShape) .
        mapPair (toChunky, map toChunky) .
        uncurry (limitMinRelativeValues (fromChunky minShape))
{-
*Synthesizer.Plain.Oscillator> let ip = Interpolation.linear in limitRelativeShapes ip ip 13 (take 100 $ iterate (1+) (0::Double)) (0::Double, cycle [0.5,1.5])
(13.0,[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5*** Exception: Numeric.NonNegative.Chunky.-: negative number
-}


limitMinRelativeValues :: (Additive.C a, Ord a) =>
   a -> a -> [a] -> (a, [a])
limitMinRelativeValues xMin x0 xs =
   let (ys,zs) =
          span ((<zero).fst) (zip (scanl (+) (x0-xMin) xs) (x0:xs))
   in  case ys of
          [] -> (x0,xs)
          (_:yr) -> (xMin, replicateMatch yr zero ++
              case zs of
                 [] -> []
                 (z:zr) -> fst z : map snd zr)

limitMaxRelativeValues :: (Additive.C a, Ord a) =>
   a -> a -> [a] -> (a, [a])
limitMaxRelativeValues xMax x0 xs =
   let (ys,zs) =
          span (>zero) (scanl (-) (xMax-x0) xs)
   in  maybe
          (xMax, replicateMatch xs zero)
          (\ ~(yl,yr) -> (x0, takeMatch yl xs ++ takeMatch zs (yr : repeat zero)))
          (viewListR ys)

{- |
Avoids negative numbers and thus can be used with Chunky numbers.
-}
limitMaxRelativeValuesNonNeg :: (Additive.C a, Ord a) =>
   a -> a -> [a] -> (a, [a])
limitMaxRelativeValuesNonNeg xMax x0 xs =
   let (ys,zs) =
          span fst (scanl (\(_,acc) d -> safeSub acc d) (safeSub xMax x0) xs)
   in  maybe
          (xMax, replicateMatch xs zero)
          (\ ~(yl, ~(_,yr)) -> (x0, takeMatch yl xs ++ takeMatch zs (yr : repeat zero)))
          (viewListR ys)
{-
*Synthesizer.Plain.Oscillator> limitMaxRelativeValuesNonNeg (let inf = 1+inf in inf) (0::Chunky.T NonNeg.Rational) (repeat 2.5)
-}

safeSub :: (Additive.C a, Ord a) => a -> a -> (Bool, a)
safeSub a b = (a>=b, a-b)