{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE Rank2Types #-}
{- |
<http://arxiv.org/abs/0911.5171>
-}
module Synthesizer.LLVM.CausalParameterized.Helix (
   -- * time and phase control based on the helix model
   static,
   staticPacked,
   dynamic,
   dynamicLimited,

   -- * useful control curves
   zigZag,
   zigZagPacked,
   zigZagLong,
   zigZagLongPacked,
   ) where

import qualified Synthesizer.LLVM.CausalParameterized.ProcessValue as CausalPV
import qualified Synthesizer.LLVM.CausalParameterized.ProcessPacked as CausalPS
import qualified Synthesizer.LLVM.CausalParameterized.ProcessPrivate
                                                              as CausalPrivP
import qualified Synthesizer.LLVM.CausalParameterized.Process as CausalP
import qualified Synthesizer.LLVM.CausalParameterized.Functional as Func
import qualified Synthesizer.LLVM.Parameterized.SignalPacked as SigPS
import qualified Synthesizer.LLVM.Parameterized.SignalPrivate as SigP
import qualified Synthesizer.LLVM.CausalParameterized.RingBufferForward as RingBuffer
import qualified Synthesizer.LLVM.Frame.SerialVector as Serial
import qualified Synthesizer.LLVM.Simple.Value as Value
import qualified Synthesizer.LLVM.Interpolation as Ip
import Synthesizer.LLVM.CausalParameterized.Functional (($&), (&|&))
import Synthesizer.LLVM.CausalParameterized.Process (($*), ($<))
import Synthesizer.LLVM.Simple.Value ((%>), (%>=), (?), (??))

import qualified Synthesizer.LLVM.Storable.Vector as SVU
import qualified Data.StorableVector as SV

import qualified LLVM.DSL.Parameter as Param

import qualified LLVM.Extra.ScalarOrVector as SoV
import qualified LLVM.Extra.Vector as Vector
import qualified LLVM.Extra.Arithmetic as A
import qualified LLVM.Extra.Storable as Storable
import qualified LLVM.Extra.Marshal as Marshal
import qualified LLVM.Extra.Memory as Memory
import qualified LLVM.Extra.Tuple as Tuple
import qualified LLVM.Extra.MaybeContinuation as MaybeCont

import qualified LLVM.Core as LLVM
import LLVM.Core (CodeGenFunction, Value, IsSized, IsFloating)

import qualified Type.Data.Num.Decimal as TypeNum

import Data.Word (Word)

import Foreign.ForeignPtr (touchForeignPtr)

import Control.Arrow (first, (<<<), (^<<), (<<^))
import Control.Category (id)
import Control.Applicative (liftA2)
import Control.Functor.HT (unzip)
import Data.Traversable (mapM)
import Data.Tuple.HT (mapFst)

import qualified Algebra.Field as Field
import qualified Algebra.Ring as Ring

import NumericPrelude.Numeric hiding (splitFraction)
import NumericPrelude.Base hiding (unzip, zip, mapM, id)

import Prelude ()


{- |
Inputs are @(shape, phase)@.

The shape parameter is limited at the beginning and at the end
such that only available data is used for interpolation.
Actually, we allow almost one step less than possible,
since the right boundary of the interval of admissible @shape@ values is open.
-}
static ::
   (Storable.C vh, Tuple.ValueOf vh ~ v,
    Ip.C nodesStep, Ip.C nodesLeap,
    SoV.RationalConstant a, SoV.Fraction a,
    Marshal.C a, Tuple.ValueOf a ~ Value a, LLVM.IsPrimitive a) =>
   (forall r. Ip.T r nodesLeap (Value a) v) ->
   (forall r. Ip.T r nodesStep (Value a) v) ->
   Param.T p Int ->
   Param.T p a ->
   Param.T p (SV.Vector vh) ->
   CausalP.T p (Value a, Value a) v
static ipLeap ipStep periodInt period vec =
   let period32 = Param.wordInt periodInt
       cellMargin = combineMarginParams ipLeap ipStep periodInt
   in  interpolateCell ipLeap ipStep
       <<<
       first (peekCell cellMargin period32 vec)
       <<<
       flattenShapePhaseProc period32 period
       <<<
       first
          (limitShape cellMargin period32
              (Param.wordInt $ fmap SV.length vec))


staticPacked ::
   (Storable.C vh, Tuple.ValueOf vh ~ ve, Serial.Element v ~ ve,
    Ip.C nodesStep, Ip.C nodesLeap,
    Serial.Size (nodesLeap (nodesStep v)) ~ n,
    Serial.C (nodesLeap (nodesStep v)),
    Serial.Element (nodesLeap (nodesStep v)) ~
       nodesLeap (nodesStep (Serial.Element v)),
    TypeNum.Positive n,
    SoV.RationalConstant a, SoV.Fraction a, Vector.Real a,
    Marshal.C a, Tuple.ValueOf a ~ Value a, LLVM.IsPrimitive a) =>
   (forall r. Ip.T r nodesLeap (Serial.Value n a) v) ->
   (forall r. Ip.T r nodesStep (Serial.Value n a) v) ->
   Param.T p Int ->
   Param.T p a ->
   Param.T p (SV.Vector vh) ->
   CausalP.T p (Serial.Value n a, Serial.Value n a) v
staticPacked ipLeap ipStep periodInt period vec =
   let period32 = Param.wordInt periodInt
       cellMargin = combineMarginParams ipLeap ipStep periodInt
   in  interpolateCell ipLeap ipStep
       <<<
       first (CausalPS.pack
          (peekCell (fmap elementMargin cellMargin) period32 vec))
       <<<
       flattenShapePhaseProcPacked period32 period
       <<<
       first
          (limitShapePacked cellMargin period32
              (Param.wordInt $ fmap SV.length vec))


{- |
In contrast to 'dynamic' this one ends
when the end of the manipulated signal is reached.
-}
dynamicLimited ::
   (Ip.C nodesStep, Ip.C nodesLeap,
    A.Additive v, Memory.C v,
    SoV.RationalConstant a, SoV.Fraction a,
    Marshal.C a, Tuple.ValueOf a ~ Value a, LLVM.IsPrimitive a, LLVM.CmpRet a) =>
   (forall r. Ip.T r nodesLeap (Value a) v) ->
   (forall r. Ip.T r nodesStep (Value a) v) ->
   Param.T p Int ->
   Param.T p a ->
   SigP.T p v ->
   CausalP.T p (Value a, Value a) v
dynamicLimited ipLeap ipStep periodInt period sig =
   dynamicGen
      (\cellMargin (skips, fracs) ->
         let windows =
                RingBuffer.trackSkip (fmap Ip.marginNumber cellMargin) sig $& skips
         in  (windows,
              CausalP.delay1Zero $& skips,
              CausalP.delay1Zero $& fracs))
      ipLeap ipStep periodInt period

{- |
If the time control exceeds the end of the input signal,
then the last waveform is locked.
This is analogous to 'static'.
-}
dynamic ::
   (Ip.C nodesStep, Ip.C nodesLeap,
    A.Additive v, Memory.C v,
    SoV.RationalConstant a, SoV.Fraction a,
    Marshal.C a, Tuple.ValueOf a ~ Value a, LLVM.IsPrimitive a, LLVM.CmpRet a) =>
   (forall r. Ip.T r nodesLeap (Value a) v) ->
   (forall r. Ip.T r nodesStep (Value a) v) ->
   Param.T p Int ->
   Param.T p a ->
   SigP.T p v ->
   CausalP.T p (Value a, Value a) v
dynamic ipLeap ipStep periodInt period sig =
   dynamicGen
      (\cellMargin (skips, fracs) ->
         let {-
             For conformance with 'static'
             we stop one step before the definite end.
             We achieve this by using a buffer
             that is one step longer than necessary.
             -}
             ((running, actualSkips), windows) =
                mapFst unzip $ unzip $
                RingBuffer.trackSkipHold
                   (fmap (succ . Ip.marginNumber) cellMargin) sig $& skips
             holdFracs =
                CausalPV.zipWithSimple (\r fr -> r ? (fr, 1))
                $&
                running &|& (CausalP.delay1Zero $& fracs)
         in  (windows, actualSkips, holdFracs))
      ipLeap ipStep periodInt period

dynamicGen ::
   (Ip.C nodesStep, Ip.C nodesLeap,
    A.Additive v, Memory.C v,
    SoV.RationalConstant a, SoV.Fraction a,
    Marshal.C a, Tuple.ValueOf a ~ Value a, LLVM.IsPrimitive a, LLVM.CmpRet a) =>
   (Param.T p (Ip.Margin (nodesLeap (nodesStep v))) ->
    (Func.T p (Value a, Value a) (Value Word),
     Func.T p (Value a, Value a) (Value a)) ->
    (Func.T p (Value a, Value a) (RingBuffer.T v),
     Func.T p (Value a, Value a) (Value Word),
     Func.T p (Value a, Value a) (Value a))) ->
   (forall r. Ip.T r nodesLeap (Value a) v) ->
   (forall r. Ip.T r nodesStep (Value a) v) ->
   Param.T p Int ->
   Param.T p a ->
   CausalP.T p (Value a, Value a) v
dynamicGen limitMaxShape ipLeap ipStep periodInt period =
   let period32 = Param.wordInt periodInt
       cellMargin = combineMarginParams ipLeap ipStep periodInt
       minShape =
          Param.wordInt $ fmap fst $
          liftA2 shapeMargin cellMargin periodInt

   in  Func.withArgs $ \(shape, phase) ->
          let (windows, skips, fracs) =
                 limitMaxShape cellMargin $
                 unzip (integrateFrac $& (limitMinShape minShape $& shape))
              (offsets, shapePhases) =
                 unzip
                    (flattenShapePhaseProc period32 period $&
                       (constantFromWord32 minShape + fracs)
                       &|&
                       (CausalP.osciCoreSync $&
                          phase
                          &|&
                          negate
                             (CausalPV.map (flip (/)) period $&
                                (CausalP.mapSimple LLVM.inttofp $& skips))))
          in  interpolateCell ipLeap ipStep $&
                 (CausalP.map (uncurry . cellFromBuffer) period32
                  $&
                  windows
                  &|&
                  offsets)
                 &|&
                 shapePhases

constantFromWord32 ::
   (IsFloating a, LLVM.IsPrimitive a) =>
   Param.T p Word -> Func.T p inp (Value a)
constantFromWord32 x =
   Func.fromSignal
      (CausalP.mapSimple LLVM.inttofp $* SigP.constant x)

limitMinShape ::
   (IsFloating a, IsSized a, LLVM.IsPrimitive a, LLVM.CmpRet a) =>
   Param.T p Word ->
   CausalP.T p (Value a) (Value a)
limitMinShape xLim =
   CausalPV.mapAccum
      (\_ x lim -> (x%>=lim) ? ((x-lim,zero), (zero,lim-x)))
      (Value.lift1 LLVM.inttofp) (return ()) xLim

integrateFrac ::
   (IsFloating a, IsSized a, LLVM.IsPrimitive a) =>
   CausalP.T p (Value a) (Value Word, Value a)
integrateFrac =
   CausalP.mapAccumSimple
      (\a (_n,frac) -> do
         s <- splitFraction =<< A.add a frac
         return (s, s))
      (return (A.zero, A.zero))


interpolateCell ::
   (Ip.C nodesStep, Ip.C nodesLeap) =>
   (forall r. Ip.T r nodesLeap a v) ->
   (forall r. Ip.T r nodesStep a v) ->
   CausalP.T p (nodesLeap (nodesStep v), (a, a)) v
interpolateCell ipLeap ipStep =
   CausalP.mapSimple
      (\(nodes, (leap,step)) ->
         ipLeap leap =<< mapM (ipStep step) nodes)

cellFromBuffer ::
   (Memory.C a, Ip.C nodesLeap, Ip.C nodesStep) =>
   Value Word ->
   RingBuffer.T a ->
   Value Word ->
   CodeGenFunction r (nodesLeap (nodesStep a))
cellFromBuffer periodInt buffer offset =
   Ip.indexNodes
      (Ip.indexNodes (flip RingBuffer.index buffer) A.one)
      periodInt offset

elementMargin ::
   Ip.Margin (nodesLeap (nodesStep v)) ->
   Ip.Margin (nodesLeap (nodesStep (Serial.Element v)))
elementMargin (Ip.Margin x y) = Ip.Margin x y

peekCell ::
   (Storable.C a, Tuple.ValueOf a ~ value, Ip.C nodesLeap, Ip.C nodesStep) =>
   Param.T p (Ip.Margin (nodesLeap (nodesStep value))) ->
   Param.T p Word ->
   Param.T p (SV.Vector a) ->
   CausalP.T p (Value Word) (nodesLeap (nodesStep value))
peekCell margin period32 vec =
   Param.withValue (Param.wordInt $ fmap Ip.marginOffset margin) $ \getOffset valueOffset ->
   Param.withValue period32 $ \getPeriod valuePeriod -> CausalPrivP.Cons
      (\(p,off,per) () n () -> MaybeCont.lift $ do
         offset <- LLVM.bitcast =<< A.sub n (valueOffset off)
         perInt <- LLVM.bitcast $ valuePeriod per
         nodes <-
            Ip.loadNodes (Ip.loadNodes Storable.load A.one) perInt
               =<< Storable.advancePtr offset p
         return (nodes, ()))
      (return ())
      (return . flip (,) ())
      (const $ const $ return ())
      (\p ->
         let (fp,ptr,_l) = SVU.unsafeToPointers $ Param.get vec p
         in  return (fp, (ptr, getOffset p, getPeriod p)))
      touchForeignPtr


flattenShapePhaseProc ::
   (IsFloating a, SoV.Fraction a, SoV.RationalConstant a,
    Marshal.C ah, Tuple.ValueOf ah ~ Value a, LLVM.IsPrimitive a) =>
   Param.T p Word ->
   Param.T p ah ->
   CausalP.T p
      (Value a, Value a)
      (Value Word, (Value a, Value a))
flattenShapePhaseProc period32 period =
   CausalP.map
      (\(perInt, per) (shape, phase) ->
         flattenShapePhase perInt per shape phase)
      (liftA2 (,) period32 period)

flattenShapePhaseProcPacked ::
   (IsFloating a, Vector.Real a, SoV.RationalConstant a,
    Marshal.C ah, Tuple.ValueOf ah ~ Value a, LLVM.IsPrimitive a,
    TypeNum.Positive n) =>
   Param.T p Word ->
   Param.T p ah ->
   CausalP.T p
      (Serial.Value n a, Serial.Value n a)
      (Serial.Value n Word,
       (Serial.Value n a, Serial.Value n a))
flattenShapePhaseProcPacked period32 period =
   CausalP.map
      (\(perInt, per) (Serial.Cons shape, Serial.Cons phase) -> do
         perIntVec <- SoV.replicate perInt
         perVec <- SoV.replicate per
         (i, (leap, step)) <-
            flattenShapePhase perIntVec perVec shape phase
         return (Serial.Cons i, (Serial.Cons leap, Serial.Cons step)))
      (liftA2 (,) period32 period)

flattenShapePhase ::
   (IsFloating a, SoV.Fraction a, SoV.RationalConstant a,
    LLVM.ShapeOf a ~ LLVM.ShapeOf i, LLVM.IsInteger i) =>
   Value i ->
   Value a ->
   Value a -> Value a ->
   CodeGenFunction r (Value i, (Value a, Value a))
flattenShapePhase = Value.unlift4 $ \periodInt period shape phase ->
   let qLeap = Value.lift1 A.fraction $ shape/period - phase
       (n,qStep) =
          unzip $ Value.lift1 splitFraction $
          {-
          If 'shape' is correctly limited,
          the value is always non-negative algebraically,
          but maybe not numerically.
          -}
          Value.max zero $
          shape - qLeap * Value.lift1 LLVM.inttofp periodInt
   in  (n,(qLeap,qStep))

{- |
You must make sure, that the argument is non-negative.
-}
splitFraction ::
   (IsFloating a, LLVM.IsInteger i, LLVM.ShapeOf a ~ LLVM.ShapeOf i) =>
   Value a -> CodeGenFunction r (Value i, Value a)
splitFraction x = do
   n <- LLVM.fptoint x
   frac <- A.sub x =<< LLVM.inttofp n
   return (n, frac)


limitShape ::
   (IsSized t, IsFloating t, SoV.Real t,
    LLVM.ShapeOf t ~ LLVM.ShapeOf i,
    Marshal.C i, Tuple.ValueOf i ~ Value i,
    Ring.C i, LLVM.IsInteger i, SoV.IntegerConstant i,
    Ip.C nodesStep, Ip.C nodesLeap) =>
   Param.T p (Ip.Margin (nodesLeap (nodesStep value))) ->
   Param.T p i ->
   Param.T p i ->
   CausalP.T p (Value t) (Value t)
limitShape margin periodInt len =
   CausalPV.zipWithSimple (Value.limit . unzip)
   $<
   limitShapeSignal margin periodInt len

limitShapePacked ::
   (IsSized t, IsFloating t, LLVM.IsPrimitive t, Vector.Real t,
    TypeNum.Positive n,
    Ip.C nodesStep, Ip.C nodesLeap) =>
   Param.T p (Ip.Margin (nodesLeap (nodesStep value))) ->
   Param.T p Word ->
   Param.T p Word ->
   CausalP.T p (Serial.Value n t) (Serial.Value n t)
limitShapePacked margin periodInt len =
   CausalPV.zipWithSimple
      (\minmax shape ->
         let (minShape,maxShape) = unzip minmax
         in  Value.limit
                (Value.lift1 Serial.upsample minShape,
                 Value.lift1 Serial.upsample maxShape)
                shape)
   $<
   limitShapeSignal margin periodInt len

limitShapeSignal ::
   (IsSized t, IsFloating t,
    LLVM.ShapeOf t ~ LLVM.ShapeOf i,
    Marshal.C i, Tuple.ValueOf i ~ Value i,
    Ring.C i, LLVM.IsInteger i, SoV.IntegerConstant i,
    Ip.C nodesStep, Ip.C nodesLeap) =>
   Param.T p (Ip.Margin (nodesLeap (nodesStep value))) ->
   Param.T p i ->
   Param.T p i ->
   SigP.T p (Value t, Value t)
limitShapeSignal margin periodInt len =
   SigP.Cons
      (\minMax () () -> return (minMax, ()))
      (return ())
      (\(minShapeInt, maxShapeInt) -> do
         minShape <- LLVM.inttofp minShapeInt
         maxShape <- LLVM.inttofp maxShapeInt
         return ((minShape, maxShape), ()))
      (const $ const $ return ())
      (\p -> return ((),
         shapeLimits
            (Param.get margin p)
            (Param.get periodInt p)
            (Param.get len p)))
      (const $ return ())


_limitShape ::
   (Ring.C th, Marshal.C th, Tuple.ValueOf th ~ t, A.Real t,
    Ip.C nodesStep, Ip.C nodesLeap) =>
   Ip.Margin (nodesLeap (nodesStep value)) ->
   Param.T p th ->
   Param.T p th ->
   CausalP.T p t t
_limitShape margin periodInt len =
   CausalPrivP.Cons
      (\(minShape,maxShape) () shape () -> MaybeCont.lift $ do
         limited <- A.min maxShape =<< A.max minShape shape
         return (limited, ()))
      (return ())
      (\minmax -> return (minmax, ()))
      (const $ const $ return ())
      (\p ->
         return
            ((),
             shapeLimits margin
                (Param.get periodInt p)
                (Param.get len p)))
      (const $ return ())

shapeLimits ::
   (Ip.C nodesLeap, Ip.C nodesStep, Ring.C t) =>
   Ip.Margin (nodesLeap (nodesStep value)) ->
   t ->
   t ->
   (t, t)
shapeLimits margin periodInt len =
   case shapeMargin margin periodInt of
      (leftMargin, rightMargin) ->
         (leftMargin, len - rightMargin)

_shapeLimits ::
   (Ip.C nodesLeap, Ip.C nodesStep,
    IsFloating t, LLVM.ShapeOf t ~ LLVM.ScalarShape) =>
   Ip.Margin (nodesLeap (nodesStep value)) ->
   Value.T (Value Word) ->
   Value.T (Value t) ->
   (Value.T (Value t), Value.T (Value t))
_shapeLimits margin periodInt len =
   let (leftMargin, rightMargin) = shapeMargin margin periodInt
   in  (Value.lift1 LLVM.inttofp leftMargin,
        len - Value.lift1 LLVM.inttofp rightMargin)

shapeMargin ::
   (Ip.C nodesLeap, Ip.C nodesStep, Ring.C i) =>
   Ip.Margin (nodesLeap (nodesStep value)) ->
   i -> (i, i)
shapeMargin margin periodInt =
   let leftMargin = fromIntegral (Ip.marginOffset margin) + periodInt
       rightMargin = fromIntegral (Ip.marginNumber margin) - leftMargin
   in  (leftMargin, rightMargin)

combineMarginParams ::
   (Ip.C nodesStep, Ip.C nodesLeap) =>
   (forall r. Ip.T r nodesLeap a v) ->
   (forall r. Ip.T r nodesStep a v) ->
   Param.T p Int ->
   Param.T p (Ip.Margin (nodesLeap (nodesStep v)))
combineMarginParams ipLeap ipStep periodInt =
   fmap
      (combineMargins (Ip.toMargin ipLeap) (Ip.toMargin ipStep))
      periodInt

combineMargins ::
   Ip.Margin (nodesLeap value) ->
   Ip.Margin (nodesStep value) ->
   Int ->
   Ip.Margin (nodesLeap (nodesStep value))
combineMargins marginLeap marginStep periodInt =
   Ip.Margin {
      Ip.marginNumber =
         Ip.marginNumber marginStep +
         Ip.marginNumber marginLeap * periodInt,
      Ip.marginOffset =
         Ip.marginOffset marginStep +
         Ip.marginOffset marginLeap * periodInt
   }


{- |
@zigZagLong loopStart loopLength@
creates a curve that starts at 0
and is linear until it reaches @loopStart+loopLength@.
Then it begins looping in a ping-pong manner
between @loopStart+loopLength@ and @loopStart@.
It is useful as @shape@ control for looping a sound.
Input of the causal process is the slope (or frequency) control.
Slope values must not be negative.

*Main> Sig.renderChunky SVL.defaultChunkSize (Causal.take 25 <<< Helix.zigZagLong 6 10 $* 2) () :: SVL.Vector Float
VectorLazy.fromChunks [Vector.pack [0.0,1.999999,3.9999995,6.0,8.0,10.0,12.0,14.0,15.999999,14.000001,12.0,10.0,7.999999,6.0,8.0,10.0,12.0,14.0,16.0,14.0,11.999999,9.999998,7.999998,6.0000024,8.000002]]
-}
zigZagLong ::
   (Marshal.C a, Tuple.ValueOf a ~ Value a,
    SoV.Fraction a, IsFloating a, SoV.RationalConstant a, LLVM.CmpRet a,
    Field.C a) =>
   Param.T p a ->
   Param.T p a ->
   CausalP.T p (Value a) (Value a)
zigZagLong =
   zigZagLongGen (CausalP.fromSignal . SigP.constant) zigZag

zigZagLongPacked ::
   (Marshal.C a, Tuple.ValueOf a ~ Value a,
    Marshal.Vector n a, Tuple.VectorValueOf n a ~ Value (LLVM.Vector n a),
    SoV.Fraction a, SoV.RationalConstant a, Vector.Real a,
    LLVM.IsPrimitive a, Field.C a,
    (n TypeNum.:*: LLVM.SizeOf a) ~ asize,
    TypeNum.Positive asize,
    TypeNum.Positive n) =>
   Param.T p a ->
   Param.T p a ->
   CausalP.T p (Serial.Value n a) (Serial.Value n a)
zigZagLongPacked =
   zigZagLongGen (CausalP.fromSignal . SigPS.constant) zigZagPacked

zigZagLongGen ::
   (A.RationalConstant al, A.Field al, Field.C a) =>
   (Param.T p a -> CausalP.T p al al) ->
   (Param.T p a -> CausalP.T p al al) ->
   Param.T p a ->
   Param.T p a ->
   CausalP.T p al al
zigZagLongGen constant zz prefix loop =
   zz (negate $ prefix/loop) * constant loop + constant prefix
   <<<
   id / constant loop

{- |
@zigZag start@ creates a zig-zag curve with values between 0 and 1, inclusively,
that is useful as @shape@ control for looping a sound.
Input of the causal process is the slope (or frequency) control.
Slope values must not be negative.
The start value must be at most 2 and may be negative.
-}
zigZag ::
   (Marshal.C a, Tuple.ValueOf a ~ Value a,
    SoV.Fraction a, IsFloating a, SoV.RationalConstant a, LLVM.CmpRet a) =>
   Param.T p a ->
   CausalP.T p (Value a) (Value a)
zigZag start =
   CausalPV.mapSimple (\x -> 1-abs (1-x))
   <<<
   CausalPV.mapAccum
      (\_ d t0 ->
         let t1 = t0+d
         in  (t0, wrap (curry . (??)) t1))
      id (return ()) start

zigZagPacked ::
   (Marshal.C a, Tuple.ValueOf a ~ Value a,
    SoV.Fraction a, IsFloating a, Vector.Real a, SoV.RationalConstant a,
    LLVM.CmpRet a,
    TypeNum.Positive n) =>
   Param.T p a ->
   CausalP.T p (Serial.Value n a) (Serial.Value n a)
zigZagPacked start =
   Serial.Cons
   ^<<
   CausalPV.mapSimple (\x -> 1 - abs (1-x))
   <<<
   CausalPV.mapAccum
      (\_ d t0 ->
         let (t1, cum) = unzip $ Value.lift2 Vector.cumulate t0 d
         {-
         LLVM.select can be replaced by (??)
         once vector select is implemented by LLVM.
         -}
         in  (wrap (Value.lift3 LLVM.select) cum, t1))
      id (return ()) start
   <<^
   (\(Serial.Cons v) -> v)

wrap ::
   (SoV.RationalConstant a, IsFloating a, SoV.Fraction a, LLVM.CmpRet a) =>
   (Value.T (Value (LLVM.CmpResult a)) ->
    Value.T (Value a) ->
    Value.T (Value a) ->
    Value.T (Value a)) ->
   Value.T (Value a) -> Value.T (Value a)
wrap select a = select (a%>0) (2 * Value.fraction (a/2)) a