-- Copyright 2017 Evan Laforge
-- This program is distributed under the terms of the GNU General Public
-- License 3.0, see COPYING or http://www.gnu.org/licenses/gpl-3.0.txt

{-# LANGUAGE CPP #-}
-- | The 'Signal' type and functions.
module Util.Segment (
    Signal, SignalS, Interpolate
    , Segment(..)
    , X, Sample(..)
    -- * construct / destruct
    , empty
    , constant, constant_val, constant_val_num
    , beginning
    , from_vector, to_vector
    , from_samples, to_samples, to_samples_desc
    , from_pairs, to_pairs, to_pairs_desc
    , from_segments, to_segments, samples_to_segments
    , simplify
    , unfoldr
    , with_ptr

    -- * query
    , null
    , at, at_negative, segment_at
    , head, last
    , maximum, minimum
    , find

    -- * concat
    , concat, prepend
    -- * slice
    , drop_after, clip_after, num_clip_after
    , drop_before, clip_before, clip_before_samples
    -- * transform
    , shift
    , map_y, map_y_linear, map_x
    -- , map_segments
    , transform_samples, map_err

    -- ** hacks
    , drop_discontinuity_at

    -- * Boxed
    , Boxed
    -- * NumSignal
    , NumSignal
    , num_interpolate, num_interpolate_s
    , invert
    , integrate
    -- * resample
    , linear_operator
    , resample_num, resample_maybe
    , sample_xs, add_zero_transition
    -- * piecewise constant
    , to_piecewise_constant
#ifdef TESTING
    , module Util.Segment
#endif
) where
import           Prelude hiding (concat, head, last, maximum, minimum, null)
import qualified Control.DeepSeq as DeepSeq
import qualified Data.List as List
import qualified Data.Maybe as Maybe
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Storable as Vector.Storable

import qualified Foreign

import qualified Util.Pretty as Pretty
import qualified Util.Seq as Seq
import qualified Util.Serialize as Serialize
import qualified Util.TimeVector as TimeVector
import           Util.TimeVector (X, Sample(..))
import qualified Util.Vector

import qualified Perform.RealTime as RealTime
import qualified Ui.Types as Types

import           Global


{- | A signal modeled as segments.  Presumably the segments are linear, but
    since you pass your own 'Interpolate', nothing in this module enforces
    that, though there are some transformations that are only valid for linear
    segments.

    A signal has no value before its first sample, and maintains a constant
    value of the last sample forever.  There is an implicit discontinuity
    to the first sample, so if @x@ is the first sample, then @[(x, y), ..]@ is
    implicitly @[(x, Nothing), (x, y), ...]@.  'NumSignal' uses 0 for not set,
    so unless the first @y@ is also 0 it becomes @[(x, 0), (x, y), ..]@.

    This comes with a built-in X offset, so translation is cheap, via 'shift'.

    Each X should be >= the previous X, and there shouldn't be more than two
    equal Xs in a row.  The first ensures that binary search works, and the
    second insures that I don't try to interpolate a zero length segment.
    Construction via 'from_samples' should establish them, and transformations
    should maintain them.

    However, a few functions in here can break them, e.g. 'map_x' and 'invert',
    and I think trying to fix them would be expensive.  So be careful with
    those.  Functions should be robust against zero length segments, but if you
    break ordering you're out of luck.

    If the @y@ value doesn't have an Eq instance, there's no way to filter
    out redundant segments like @[(0, 1), (1, 1), (2, 1)]@.  Functions
    specialized to 'NumSignal' may make some effort to do that, but only if it
    seems worth it.
-}
data Signal v = Signal {
    _offset :: !X
    , _vector :: !v
    } deriving (Eq, Show)

type SignalS v y = Signal (v (Sample y))

type Interpolate y = Sample y -> Sample y -> X -> y

instance Pretty v => Pretty (Signal v) where
    format (Signal offset vector) = "Signal" Pretty.<+> Pretty.text offset_s
        Pretty.<+> Pretty.format vector
        where
        offset_s
            | offset > 0 = "+" <> pretty offset
            | otherwise = pretty offset

instance Serialize.Serialize v => Serialize.Serialize (Signal v) where
    put (Signal offset vec) = Serialize.put offset >> Serialize.put vec
    get = Signal <$> Serialize.get <*> Serialize.get

instance DeepSeq.NFData v => DeepSeq.NFData (Signal v) where
    rnf (Signal offset vec) = DeepSeq.rnf offset `seq` DeepSeq.rnf vec `seq` ()

modify_vector :: (a -> b) -> Signal a -> Signal b
modify_vector modify sig = sig { _vector = modify (_vector sig) }

data Segment y = Segment {
    _x1 :: !X, _y1 :: !y
    , _x2 :: !X, _y2 :: !y
    } deriving (Eq, Show)

instance Pretty y => Pretty (Segment y) where
    pretty (Segment x1 y1 x2 y2) = "(" <> pretty x1 <> ", " <> pretty y1
        <> ")--(" <> pretty x2 <> ", " <> pretty y2 <> ")"

-- * construct / destruct

empty :: V.Vector v a => Signal (v a)
empty = Signal 0 V.empty

constant :: V.Vector v (Sample y) => y -> SignalS v y
constant a = from_vector $ V.fromList [Sample beginning a]

constant_val :: V.Vector v (Sample a) => SignalS v a -> Maybe a
constant_val sig = case TimeVector.uncons (_vector sig) of
    -- This will naturally disregard 'shift's, which is as it should be for
    -- so-called constant signals.
    Just (Sample x1 y1, rest) | x1 <= -RealTime.large && V.null rest ->
        Just y1
    _ -> Nothing

-- | 'constant_val' for 'NumSignal's can be more clever, because it can compare
-- Ys.  Also NumSignals are implicitly 0 before the first sample.
constant_val_num :: X -> NumSignal -> Maybe Y
constant_val_num from sig = case TimeVector.uncons (_vector sig) of
    -- I compare multiple samples because a track might have redundant
    -- values, but I still want to detect if it's constant.
    Just (Sample x y, rest)
        | x <= (from - _offset sig) && V.all ((==y) . sy) rest -> Just y
        | V.all ((==0) . sy) (_vector sig) -> Just 0
        | otherwise -> Nothing
    Nothing -> Just 0

-- | Use this as the stand-in for "since the beginning of time."
beginning :: RealTime.RealTime
beginning = -RealTime.larger

from_vector :: v -> Signal v
from_vector = Signal 0

to_vector :: V.Vector v (Sample y) => SignalS v y -> v (Sample y)
to_vector sig
    | offset == 0 = _vector sig
    | otherwise = TimeVector.map_x (+offset) (_vector sig)
    where offset = _offset sig

-- | The final sample extends for "all time".  However, there's no value before
-- the first sample.  The reason is that I'd have to have a zero value for y,
-- and there isn't really an appropriate one for pitch.
--
-- TODO I could simplify straight lines, but then I'd need Eq on y.  Maybe do
-- that separately for NumSignal.
from_samples :: V.Vector v (Sample y) => [Sample y] -> SignalS v y
from_samples =
    from_vector . V.fromList . drop_coincident
        . drop_initial_dup
        . Maybe.catMaybes . snd . List.mapAccumL in_order Nothing
    where
    -- Since the first sample comes from 0, I can drop leading dups.
    drop_initial_dup (s1 : ss@(s2 : _)) | sx s1 == sx s2 = drop_initial_dup ss
    drop_initial_dup s = s
    -- Drop out-of-order samples.
    in_order Nothing cur = (Just cur, Just cur)
    in_order (Just prev) cur
        | sx cur < sx prev = (Just prev, Nothing)
        | otherwise = (Just cur, Just cur)
    -- Abbreviate coincident samples.
    drop_coincident (Sample x1 y1 : _ : sn@(Sample x2 _ : _)) | x1 == x2 =
        drop_coincident $ Sample x1 y1 : sn
    drop_coincident (s1:sn) = s1 : drop_coincident sn
    drop_coincident [] = []

to_samples :: V.Vector v (Sample y) => SignalS v y -> [Sample y]
to_samples sig =
    (if _offset sig == 0 then id else map (plus (_offset sig))) $
        V.toList (_vector sig)
    where
    plus n (Sample x y) = Sample (n+x) y
-- to_samples = V.toList . to_vector
-- TODO verify that TimeVector.map_x fuses with V.toList so there is no extra
-- vector.

to_samples_desc :: V.Vector v (Sample y) => SignalS v y -> [Sample y]
to_samples_desc (Signal offset vec) =
    (if offset == 0 then id else map (plus offset)) $
        Util.Vector.to_reverse_list vec
    where
    plus n (Sample x y) = Sample (n+x) y

from_pairs :: V.Vector v (Sample y) => [(X, y)] -> SignalS v y
from_pairs = from_samples . map (uncurry Sample)

to_pairs :: V.Vector v (Sample y) => SignalS v y -> [(X, y)]
to_pairs = map TimeVector.to_pair . to_samples

to_pairs_desc :: V.Vector v (Sample y) => SignalS v y -> [(X, y)]
to_pairs_desc = map TimeVector.to_pair . to_samples_desc

from_segments :: V.Vector v (Sample y) => [Segment y] -> SignalS v y
from_segments = from_samples . to_list
    where
    to_list (Segment x1 y1 x2 y2 : segments) =
        Sample x1 y1 : Sample x2 y2 : to_list segments
    to_list [] = []

to_segments :: V.Vector v (Sample y) => SignalS v y -> [Segment y]
to_segments = samples_to_segments . to_samples

samples_to_segments :: [Sample y] -> [Segment y]
samples_to_segments = go
    where
    go [] = []
    go [Sample x y]
        | x < RealTime.large = [Segment x y RealTime.large y]
        | otherwise = []
    go (Sample x1 y1 : xs@(Sample x2 y2 : _))
        | x1 == x2 = go xs
        | otherwise = Segment x1 y1 x2 y2 : go xs

-- | Simplify away redundant samples.
simplify :: (Eq x, Eq y) => [(x, y)] -> [(x, y)]
simplify ((x1, _) : xys@((x2, _) : _)) | x1 == x2 = simplify xys
simplify xys = go xys
    where
    -- Drop samples in the middle of horizontal and vertical lines.
    go ((x1, y1) : (x2, y2) : xys@((x3, y3) : _))
        | y1 == y2 && y2 == y3 = go ((x1, y1) : xys)
        | x1 == x2 && x2 == x3 = go ((x1, y1) : xys)
    -- Identical samples are always redundant.
    go ((x1, y1) : xys@((x2, y2) : _))
        | x1 == x2 && y1 == y2 = go xys
    go xys = xys
{-# INLINEABLE simplify #-}
{-# SPECIALIZE simplify :: [(X, Y)] -> [(X, Y)] #-}

unfoldr :: V.Vector v (Sample y) => (state -> Maybe ((X, y), state)) -> state
    -> SignalS v y
unfoldr gen state = from_vector $ TimeVector.unfoldr gen state

-- | Get a Ptr to the vector.  This is 'Vector.Storable.unsafeWith'.
with_ptr :: Foreign.Storable a =>
    Signal (Vector.Storable.Vector a) -> (X -> Foreign.Ptr a -> Int-> IO b)
    -> IO b
with_ptr (Signal offset vec) action = TimeVector.with_ptr vec (action offset)

-- * query

null :: V.Vector v (Sample y) => SignalS v y -> Bool
null = V.null . _vector

-- | The arguments may seem backwards, but I've always done it this way, and it
-- seems to be more convenient in practice.
at :: V.Vector v (Sample y) => Interpolate y -> X -> SignalS v y -> Maybe y
at interpolate x_ (Signal offset vec)
    | i < 0 = Nothing
    | i + 1 == V.length vec = Just (sy (V.unsafeIndex vec i))
    | otherwise =
        Just $ interpolate (V.unsafeIndex vec i) (V.unsafeIndex vec (i+1)) x
    where
    i = TimeVector.highest_index x vec
    x = x_ - offset

-- | Like 'at', but if the x matches a discontinuity, take the value before
-- instead of after.
at_negative :: V.Vector v (Sample y) => Interpolate y -> X -> SignalS v y
    -> Maybe y
at_negative interpolate x signal = do
    Segment x1 y1 x2 y2 <- segment_at_orientation Types.Negative x signal
    return $ interpolate (Sample x1 y1) (Sample x2 y2) x

segment_at :: V.Vector v (Sample y) => X -> SignalS v y -> Maybe (Segment y)
segment_at  = segment_at_orientation Types.Positive

segment_at_orientation :: V.Vector v (Sample y) => Types.Orientation -> X
    -> SignalS v y -> Maybe (Segment y)
segment_at_orientation orient x (Signal offset vec) =
    bump . snd <$> segment_at_v orient (x - offset) vec
    where bump (Segment x1 y1 x2 y2) = Segment (x1+offset) y1 (x2+offset) y2

segment_at_v :: V.Vector v (Sample y) => Types.Orientation -> X -> v (Sample y)
    -> Maybe (Int, Segment y)
segment_at_v orient x vec
    | i < 0 = Nothing
    | otherwise =
        let Sample x1 y1 = V.unsafeIndex vec i
            Sample x2 y2 = if i + 1 >= V.length vec
                then Sample RealTime.large y1
                else V.unsafeIndex vec (i+1)
        in Just (i, Segment x1 y1 x2 y2)
    where
    i = get x vec
    get = case orient of
        Types.Negative -> TimeVector.index_below
        Types.Positive -> TimeVector.highest_index

head :: V.Vector v (Sample y) => SignalS v y -> Maybe (X, y)
head sig = case TimeVector.head (_vector sig) of
    Nothing -> Nothing
    Just (Sample x y) -> Just (_offset sig + x, y)

last :: V.Vector v (Sample y) => SignalS v y -> Maybe (X, y)
last sig = case TimeVector.last (_vector sig) of
    Nothing -> Nothing
    Just (Sample x y) -> Just (_offset sig + x, y)

minimum, maximum :: (V.Vector v (Sample a), Ord a) => SignalS v a -> Maybe a
minimum sig
    | null sig = Nothing
    | otherwise = Just $ sy $ V.minimumBy (\a b -> compare (sy a) (sy b)) $
        _vector sig
maximum sig
    | null sig = Nothing
    | otherwise = Just $ sy $ V.maximumBy (\a b -> compare (sy a) (sy b)) $
        _vector sig

find :: V.Vector v (Sample y) => (X -> y -> Bool) -> Signal (v (Sample y))
    -> Maybe (X, y)
find f sig = first (+ _offset sig) . TimeVector.to_pair <$>
    V.find (\(Sample x y) -> f (x + _offset sig) y) (_vector sig)

-- * concat

-- | Concatenate signals, where signals to the right replace the ones to the
-- left where they overlap.
concat :: V.Vector v (Sample y) => Maybe (y -> y -> Bool)
    -- ^ signals with Eq y can drop some redundant samples
    -> Interpolate y -> [SignalS v y] -> SignalS v y
concat _ _ [] = empty
concat _ _ [sig] = sig
concat maybe_eq interpolate sigs =
    from_vector . V.concat . try_strip_duplicates . reverse . chunks . reverse
        . map to_vector $ sigs
    where
    chunks [] = []
    chunks [v] = [v]
    -- head of v1 cuts of tail of v2
    -- v1:     |--->        |--->
    -- v2:   |--->        |->
    -- vs: |--->     => |->
    chunks (v1:v2:vs) = case sx <$> TimeVector.head v1 of
        Nothing -> chunks (v2:vs)
        Just x1 -> case TimeVector.last clipped of
            Nothing -> chunks (v1:vs)
            Just end
                | sx end < x1 -> v1 : extension end : chunks (clipped:vs)
                | otherwise -> v1 : chunks (clipped:vs)
            where
            clipped = clip_after_v interpolate x1 v2
            extension end = V.singleton (Sample x1 (sy end))
    try_strip_duplicates = case maybe_eq of
        Nothing -> id
        Just eq -> strip_duplicates eq
    -- If I have Eq, I can strip redundant Y values.
    strip_duplicates eq (v1 : v2 : vs)
        | Just (Sample x1 y1) <- TimeVector.last v1
        , Just (Sample x2 y2) <- TimeVector.head v2
        , x1 == x2 && eq y1 y2 =
            v1 : strip_duplicates eq (V.drop 1 v2 : vs)
    strip_duplicates eq (v1 : vs) = v1 : strip_duplicates eq vs
    strip_duplicates _ [] = []

-- | With 'concat', each signal start clips the signal to its left.  This is
-- the other way around, the final sample in the first signal is taken as its
-- end, and it replaces the start of the second signal.
prepend :: V.Vector v (Sample y) => Maybe (y -> y -> Bool) -> Interpolate y
    -> SignalS v y -> SignalS v y -> SignalS v y
prepend eq interpolate sig1 sig2 = case last sig1 of
    Nothing -> sig2
    Just (x, _) -> concat eq interpolate [sig1, clip_before interpolate x sig2]

-- * slice

-- | Drop the segments after the given time.  The last segment may overlap it.
drop_after :: V.Vector v (Sample y) => X -> SignalS v y -> SignalS v y
drop_after x sig
    | V.null v = empty
    | otherwise = Signal { _offset = _offset sig, _vector = v }
    where v = drop_after_v (x - _offset sig) (_vector sig)

drop_after_v :: V.Vector v (Sample y) => X -> v (Sample y) -> v (Sample y)
drop_after_v x vec = case vec V.!? i of
    Nothing -> V.empty
    Just (Sample x1 _) -> V.take (if x1 >= x then i+1 else i+2) vec
    where i = TimeVector.index_below x vec

-- | This is like 'drop_after', but meant to clip the signal directly on x,
-- rather than at the first sample >=x.  This means I might have to insert a
-- new sample, which means copying the signal.  This is intended to be a "drop
-- at and after", but since signals extend infinitely to the right, I can only
-- go up to x.  TODO maybe signals should go to Nothing >= the last sample?
--
-- If the signal has only a point exactly at x, then return the empty signal.
-- This is because the first sample is effectively a transition from Nothing,
-- or 0.
clip_after :: V.Vector v (Sample y) => Interpolate y -> X -> SignalS v y
    -> SignalS v y
clip_after interpolate x sig
    | V.null v = empty
    | otherwise = Signal { _offset = _offset sig, _vector = v }
    where v = clip_after_v interpolate (x - _offset sig) (_vector sig)

clip_after_v :: V.Vector v (Sample y) => Interpolate y -> X
    -> v (Sample y) -> v (Sample y)
clip_after_v interpolate x vec
    | [Sample x0 _] <- V.toList clipped, x0 == x = V.empty
    | otherwise = case TimeVector.last clipped of
        Nothing -> V.empty
        Just (Sample x2 _)
            | x < x2, Just y <- at interpolate x (from_vector vec) ->
                V.snoc (V.take (V.length clipped - 1) clipped) (Sample x y)
            | otherwise -> clipped
    where clipped = drop_after_v x vec

num_clip_after :: Bool -> X -> NumSignal -> NumSignal
num_clip_after keep_last x sig
    | V.null clipped = empty
    | [Sample x0 _] <- V.toList clipped, x0 == x = empty
    | otherwise = Signal { _offset = _offset sig, _vector = clipped }
    where clipped = num_clip_after_v keep_last (x - _offset sig) (_vector sig)

-- | 'clip_after' specialized for 'Y'.  Since it has Eq, it can do an
-- additional optimization.
num_clip_after_v :: Bool -- ^ if False, inhibit the optimization that omits
    -- the end sample if it's a flat line
    -> X -> TimeVector.Unboxed -> TimeVector.Unboxed
num_clip_after_v keep_last x vec = case segment_at_v Types.Negative x vec of
    Nothing -> vec
    Just (i, Segment x1 y1 x2 y2)
        | not keep_last && y1 == y2 -> prefix
        | otherwise -> V.snoc prefix (Sample x (TimeVector.y_at x1 y1 x2 y2 x))
        where prefix = V.take (i+1) vec

-- | Drop the segments before the given time.  The first segment will start at
-- or before the given time.
drop_before :: V.Vector v (Sample y) => X -> SignalS v y -> SignalS v y
drop_before x sig
    | V.null clipped = empty
    | otherwise = Signal { _offset = _offset sig, _vector = clipped }
    where clipped = TimeVector.drop_before (x - _offset sig) (_vector sig)

-- | Like 'drop_before', but ensure that the signal starts exactly at the given
-- time by splitting a segment that crosses it.
clip_before :: V.Vector v (Sample y) => Interpolate y -> X -> SignalS v y
    -> SignalS v y
clip_before interpolate x sig = case head clipped of
    Nothing -> empty
    Just (x1, _)
        | x1 < x, Just y <- at interpolate x sig ->
            from_vector $ V.cons (Sample x y) (V.drop 1 (to_vector clipped))
        | otherwise -> clipped
    where clipped = drop_before x sig

-- TODO is this the same as 'to_samples . clip_before'?
clip_before_samples :: V.Vector v (Sample y) => Interpolate y -> X
    -> SignalS v y -> [Sample y]
clip_before_samples interpolate x sig = case head clipped of
    Nothing -> []
    Just (x1, _)
        | x1 < x, Just y <- at interpolate x sig ->
            Sample x y : to_samples (modify_vector (V.drop 1) clipped)
        | otherwise -> to_samples clipped

    where clipped = drop_before x sig

-- * transform

-- | Shift the signal in time.
shift :: X -> Signal v -> Signal v
shift offset sig = sig { _offset = _offset sig + offset }

-- | Apply the _offset, and set it to 0.  Just for tests.
_flatten_shift :: V.Vector v (Sample y) => SignalS v y -> SignalS v y
_flatten_shift = from_vector . to_vector

-- | Map Ys.  This resamples the signal, so it's valid for a nonlinear
-- function.
map_y :: X -> (Y -> Y) -> NumSignal -> NumSignal
map_y srate f =
    from_vector . TimeVector.map_y f . to_vector . resample_rate srate

-- | Map Ys.  Only valid if the function is linear.
map_y_linear :: V.Vector v (Sample y) => (y -> y) -> SignalS v y -> SignalS v y
map_y_linear f = modify_vector $ TimeVector.map_y f

-- | Map Xs.  The slopes will definitely change unless the function is adding
-- a constant, but presumably that's what you want.
--
-- TODO this can break 'Signal' invariants.
map_x :: V.Vector v (Sample y) => (X -> X) -> SignalS v y -> SignalS v y
map_x f = modify_vector $ TimeVector.map_x f

transform_samples :: V.Vector v (Sample y) => ([Sample y] -> [Sample y])
    -> SignalS v y -> SignalS v y
transform_samples f = from_samples . f . to_samples

map_err :: V.Vector v (Sample y) => (Sample y -> Either err (Sample y))
    -> SignalS v y -> (SignalS v y, [err])
map_err f = first from_vector . TimeVector.map_err f . to_vector

-- ** hacks

-- | Drop a x1==x2 discontinuity at the given time, if there is one.
-- Used for Block.trim_controls, which is a terrible hack that I'm trying to
-- get rid of.
drop_discontinuity_at :: V.Vector v (Sample y) => X -> SignalS v y
    -> SignalS v y
drop_discontinuity_at x sig = case V.toList clipped of
    Sample x1 _ : Sample x2 _ : _ | x == x1 && x1 == x2 ->
        from_vector $ V.concat
            [ pre
            -- Insert an extra sample to avoid changing the slope.
            , case (TimeVector.last pre, TimeVector.head post) of
                (Just (Sample _ y), Just (Sample x _)) ->
                    V.singleton (Sample x y)
                _ -> V.empty
            , drop1 post
            ]
            where
            pre = TimeVector.drop_at_after x vector
            post = TimeVector.drop_before_at x vector
    _ -> sig
    where
    vector = to_vector sig
    clipped = TimeVector.drop_before_strict (x - _offset sig) (_vector sig)
    -- Drop an extra x to avoid >2 samples in the same spot.
    drop1 v = case V.toList v of
        Sample x1 _ : Sample x2 _ : _ | x1 == x2 -> V.drop 1 v
        _ -> v

-- * Boxed

type Boxed y = Signal (TimeVector.Boxed y)

-- * NumSignal

type NumSignal = Signal TimeVector.Unboxed
type Y = TimeVector.UnboxedY

num_interpolate :: Interpolate Y
num_interpolate (Sample x1 y1) (Sample x2 y2) = TimeVector.y_at x1 y1 x2 y2

num_interpolate_s :: Segment Y -> X -> Y
num_interpolate_s (Segment x1 y1 x2 y2) = TimeVector.y_at x1 y1 x2 y2

-- | Swap X and Y.  Y must be non-decreasing or this will break 'Signal'
-- invariants.
invert :: NumSignal -> NumSignal
invert sig = Signal 0 (V.map swap (_vector sig))
    where
    swap (Sample x y) =
        Sample (RealTime.seconds y) (RealTime.to_seconds (x + _offset sig))

-- | Integrate the signal.
--
-- Since the output will have more samples than the input, this needs
-- a sampling rate.  The sampling rate determines the resolution of the tempo
-- track.  So it can probably be fairly low resolution before having
-- a noticeable impact.
integrate :: X -> NumSignal -> NumSignal
integrate srate_x =
    from_samples . List.concat . snd
        . List.mapAccumL segment 0 . map to_double . to_segments
    where
    -- Integral of nx + k = nx^2 / 2 + kx
    to_double (Segment x1 y1 x2 y2) = (s x1, y1, s x2, y2)
        where s = RealTime.to_seconds
    to_sample x y = Sample (RealTime.seconds x) y
    segment accum (x1, y1, x2, y2) =
        ( f (x2 - x1)
        , if y1 == y2
            then [to_sample x1 (f 0), to_sample x2 (f (x2 - x1))]
            else [to_sample x (f (x - x1)) | x <- Seq.range' x1 x2 (1/srate)]
        )
        where
        f x = n * x^2 / 2 + y1*x + accum
        n = (y2 - y1) / (x2 - x1)
    srate = RealTime.to_seconds srate_x

-- * resample

-- | Combine two vectors with the given function.  The signals are resampled
-- to have coincident samples, assuming linear interpolation.  This only works
-- for linear functions, so the result can also be represented with linear
-- segments.
linear_operator :: (Y -> Y -> Y) -> NumSignal -> NumSignal -> NumSignal
linear_operator merge asig bsig =
    from_samples $ zipWith3 make (get_xs ()) as2 bs2
    where
    make x ay by = Sample x (merge ay by)
    as2 = resample_num (get_xs ()) as
    bs2 = resample_num (get_xs ()) bs
    (as, bs) = to_samples2 0 asig bsig
    -- The () is to prevent memoization, which should hopefully allow the
    -- intermediate list to fuse away.  Or maybe I should try to use vectors
    -- instead of lists?
    -- TODO profile
    get_xs () = sample_xs2 (map sx as) (map sx bs)

resample_num :: [X] -> [Sample Y] -> [Y]
resample_num = resample 0 num_interpolate

-- | Like 'to_samples', except the signal that starts later gets an extra
-- sample to transition from zero.
to_samples2 :: V.Vector v (Sample y) => y -> SignalS v y -> SignalS v y
    -> ([Sample y], [Sample y])
to_samples2 zero asig bsig = case (to_samples asig, to_samples bsig) of
    (as@(Sample ax _ : _), bs@(Sample bx _ : _))
        | ax < bx -> (as, Sample bx zero : bs)
        | bx < ax -> (Sample ax zero : as, bs)
    (as, bs) -> (as, bs)

-- | The output has the union of the Xs in the inputs, except where they match
-- exactly.  Discontinuities should get two Xs.
sample_xs2 :: [X] -> [X] -> [X]
sample_xs2 = go
    where
    go [] bs = bs
    go as [] = as
    go (a:as) (b:bs)
        | a == b = a : go as bs
        | a < b = a : go as (b:bs)
        | otherwise = b : go (a:as) bs

-- ** polymorphic implementation

-- | This should be the same as 'linear_operator', except using the
-- variable length functions.  I could replace linear_operator with this, but
-- I worry that it's less efficient.
_linear_operator2 :: ([Y] -> Y) -> NumSignal -> NumSignal -> NumSignal
_linear_operator2 merge asig bsig =
    -- Seq.rotate zips up the samples from each signal.
    from_samples $ zipWith make xs $ Seq.rotate $ map (resample_num xs) samples
    where
    make x ys = Sample x (merge ys)
    xs = sample_xs $ map (map sx) samples
    samples = map (add_zero_transition 0 . to_samples) [asig, bsig]

resample :: y -> Interpolate y -> [X] -> [Sample y] -> [y]
resample = resample_ id

-- | This is the same as 'resample', only for ys without a zero.
resample_maybe :: Interpolate y -> [X] -> [Sample y] -> [Maybe y]
resample_maybe = resample_ Just Nothing

{-# INLINE resample_ #-}
resample_ :: (y1 -> y2) -> y2 -> Interpolate y1 -> [X] -> [Sample y1] -> [y2]
resample_ present absent interpolate xs samples =
    snd $ List.mapAccumL get samples xs
    where
    get ss@(Sample x1 y1 : s2s@(Sample x2 y2 : _)) x
        -- If it's a discontinuity, I want to consume the sample, or I won't
        -- see the after Y.  Each discontinuity should have 2 Xs, one for
        -- before and one for after.  This is brittle and depends on
        -- 'sample_xs2' emitting two Xs for a discontinuity and 'to_samples2'
        -- adding a "from zero" discontinuity to the first sample.  But
        -- otherwise I'd have to recognize a discontinuity here and emit one,
        -- which means this would have to be concatMap, which seems
        -- inefficient.  Of course maybe the whole thing is already so
        -- inefficient it doesn't matter.
        | x == x1 = if x1 == x2 then (s2s, present y1) else (ss, present y1)
        | x >= x2 = get s2s x
        | x > x1 = (ss, present $ interpolate (Sample x1 y1) (Sample x2 y2) x)
        | otherwise = (ss, absent)
    get ss@[Sample x1 y1] x
        | x >= x1 = (ss, present y1)
        | otherwise = (ss, absent)
    get [] _ = ([], absent)

add_zero_transition :: y -> [Sample y] -> [Sample y]
add_zero_transition zero ss@(Sample x _ : _) = Sample x zero : ss
add_zero_transition _ [] = []

-- | The output has the union of the Xs in the inputs, except where they match
-- exactly.  Discontinuities should get two Xs.  This is the list version of
-- 'sample_xs2'.
sample_xs :: [[X]] -> [X]
sample_xs = go
    where
    go [] = []
    go xss_ = case Seq.minimum xs of
        Nothing -> go (map tail xss)
        Just x -> x : go (map (drop1 (==x)) xss)
        where
        xs = map List.head xss
        xss = filter (not . List.null) xss_
    drop1 f (x:xs) | f x = xs
    drop1 _ xs = xs


-- ** constant rate resamples

-- | This is like 'to_piecewise_constant', except it retains discontinuities,
-- which is important since it's used for 'map_y', which is still operating on
-- linear segments.  Or it's like 'resample', except it uses a constant rate
-- instead of [X].
resample_rate :: X -> NumSignal -> NumSignal
resample_rate srate =
    from_samples . concatMap resample . Seq.zip_next . to_samples
    where
    resample (Sample x1 y1, Nothing) = [Sample x1 y1]
    resample (Sample x1 y1, Just (Sample x2 y2))
        | y1 == y2 || x1 == x2 = [Sample x1 y1]
        | otherwise =
            [ Sample x (TimeVector.y_at x1 y1 x2 y2 x)
            | x <- Seq.range' x1 x2 (1/srate)
            ]

-- TODO possible vector implementation that might fuse.  But this
-- requires Storable (a, b).
-- resample_rate :: X -> NumSignal -> NumSignal
-- resample_rate srate =
--     from_vector . V.concatMap resample . zip_next . to_vector
--     where
--     zip_next xs = V.zip xs (V.drop 1 xs)
--     resample (Sample x1 y1, Sample x2 y2) = V.fromList
--         [Sample x1 y1, Sample x2 y2]

to_piecewise_constant :: X -> NumSignal -> TimeVector.Unboxed
to_piecewise_constant srate =
    V.fromList . Seq.drop_dups sy . Seq.drop_initial_dups sx . List.concat
        . List.unfoldr make . to_samples
    where
    make [] = Nothing
    make [Sample x y] = Just ([Sample x y], [])
    make (Sample x1 y1 : s2s@(Sample x2 y2 : _))
        | y1 == y2 = Just ([Sample x1 y1], s2s)
        | x1 >= x2 = make s2s
        | otherwise = Just (segment x1 y1 x2 y2, s2s)
    segment x1 y1 x2 y2 =
        [ Sample x (TimeVector.y_at x1 y1 x2 y2 x)
        | x <- Seq.range' x1 x2 (1/srate)
        ]