```{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}

module Numeric.LazySplines where
import Data.List
import Control.Arrow

-- | Something that can be sampled.
class Sampleable a where
at :: a -> Double -> Double
type Poly = [Double]

-- Horner's scheme for polynomial evaluation
instance Sampleable Poly where
at x v = foldr (\c val -> c + v * val) 0 x

instance Sampleable [Poly] where
at x v = poly `at` frac
where
poly = x !! int
(int,frac) = properFraction v

type PolySegment = (Double, Poly)

type Spline = [PolySegment]

type SplinePredicate =
Double -> Double -> Poly -> Double

duration :: Spline -> Double
duration = sum . map fst

maxDuration = 10000  -- an arbitrary limit

liftS :: Double -> Spline
liftS x = [(maxDuration,[x])]
instance Sampleable PolySegment where
at (_,poly) pt = poly `at` pt

instance Sampleable Spline where
-- assume pt >= 0
at spline pt = go spline pt where
go ((dur,poly):xs) v
| v <= dur  = poly `at` v
| otherwise = go xs (v - dur)
go [] _ = error \$
"Sampling spline of out bounds " ++
show spline ++ " at: " ++ show pt

deriveSpline :: Spline -> Spline
deriveSpline = map (second diff)

integrateSpline :: Spline -> Spline
integrateSpline =
snd . mapAccumL go 0 . map (second integ)
where
go :: Double -> PolySegment -> (Double, PolySegment)
go acc (dur,poly) = (v, seg)
where v = acc + poly `at` dur
seg = (dur, realToFrac acc + poly)

inSpline2 :: (Poly -> Poly -> Poly) ->
Spline -> Spline -> Spline
inSpline2 op ((xd,x):xs) ((yd,y):ys)
| xd == yd  = (xd, v) : inSpline2 op xs ys
| xd <  yd  = (xd, v) : inSpline2 op xs (y':ys)
| otherwise = (yd, v) : inSpline2 op (x':xs) ys
where
v = x `op` y
x' = splitPoly yd (xd,x)
y' = splitPoly xd (yd,y)
splitPoly d (dur,poly) = (dur - d, shiftBy d poly)
inSpline2 _ _ _ = []

shiftBy :: Double -> Poly -> Poly
shiftBy d poly = poly # [d,1]
instance Num Spline where
fromInteger = liftS . fromInteger
negate = map (second negate)
(+)    = inSpline2 (+)
(*)    = inSpline2 (*)

mapSpline :: Bool ->
(Double -> Double -> Poly -> Poly) ->
Spline ->
Spline
mapSpline _          _ []         = []
mapSpline matchFirst f (seg:segs) =
-- leave first segment unchanged
seg :
(snd \$ mapAccumL go (dur0, seg `at` dur0) segs)
where
dur0 = fst seg

go (totalDur, lastVal) (dur, poly) =
((totalDur + dur, fnc' `at` dur), fnc')
where
fnc' = (dur, poly')

poly'
-- modify translated segment
| matchFirst = f totalDur dur \$
match lastVal poly
-- translate modified segment
| otherwise  = matchScale lastVal dur \$
f totalDur dur poly

-- replace 0-degree coefficient of a Poly
match lastVal (_:xs) = lastVal:xs
match lastVal x      = [lastVal]

-- Alter the first point, preserve the point at dur
matchScale v dur poly@(x:xs) = v : map (* scale) xs
where height = poly `at` (dur - x)
diff   = x - v
scale = height / (height - diff)
matchScale v _ x = [v]

infixl 1 `trimmingTo`
trimmingTo :: Spline -> Int -> Spline
trimmingTo spline power =
mapSpline True go spline
where
go _ _ s = take power s

infixl 1 `extrapForward`
extrapForward :: Spline -> Double -> Spline
extrapForward spline delta =
mapSpline False go spline
where
go _ _ s = shiftBy delta s

scaleRest :: Poly -> Double -> Poly
scaleRest (x:xs) c = x : map (* c) xs

infixl 1 `satisfying`
satisfying :: Spline ->
(Double, SplinePredicate) ->
Spline
satisfying spline (tol, p) =
mapSpline True go spline
where
go t d fnc = findValue tol (p t d) \$
scaleRest fnc

infixl 1 `splitWhen`
splitWhen :: Spline ->
(Double, Double, SplinePredicate) ->
Spline
splitWhen spline (tol, minsize, p) =
go 0 spline
where
go _ [] = []
go t (f@(dur,poly):fs)
| doSplit   = go t (f' : f'' : fs)
| otherwise = f : go t' fs
where

-- split if predicate is not satisfied
-- and duration is not below minsize
doSplit = dur > minsize &&
abs (p t dur poly) > tol

-- two segments, each with half
-- the duration of the original
dur' = dur / 2
f'   = (dur', poly)
f''  = (dur', shiftBy dur' poly)

t' = t + dur

infixl 1 `extendWhen`
extendWhen :: Spline ->
(Double, Double, SplinePredicate) ->
Spline
extendWhen spline (tol,maxlen,p) =
go (spline `at` 0, 0) 0 spline
where
go _ _ [] = []
go (lastVal, time) chop ((oldDur,oldPoly):fs)

-- This segment has been subsumed
-- by previous extensions,
-- disregard it.
| dur <= 0  =
go (lastVal, time) (negate dur) fs

-- Attempt to extend this segment
| otherwise =
(dur', poly') :
go (lastVal', time') chop' fs
where

-- The segment as chopped
-- to reflect prior extensions
poly = shiftBy chop oldPoly
dur  = oldDur - chop

chkPred d =
(abs \$ p time d poly) < tol

-- greatest permissible duration
-- that satisfies the predicate
dur' = lastDef dur .
takeWhile chkPred .
-- list of possible durations
takeWhile (<= maxlen) \$
iterate (* 2) dur

chop' = dur' - dur
time' = time + dur'
lastVal' = poly `at` dur'
poly' = matchScale lastVal dur' poly

lastDef def []       = def
lastDef _ [x]        = x
lastDef def (x : xs) = lastDef def xs

infixl 1 `trimSmart`
trimSmart :: Spline ->
SplinePredicate ->
Spline
trimSmart spline p =
mapSpline True go spline
where
go t dur poly =
map fst .
dropWhile chkPred \$
-- polynomials
zip polys (drop 1 polys)
where
-- possible trimmed polynomials
polys = drop 2 . inits \$  poly

chkPred (x,x') = p t dur x >
p t dur x'

-- Adopted from M. Douglas McIlroy., Power series, power serious.
-- http://www.cs.dartmouth.edu/~doug/powser.html

infixr 9 #
instance Num Poly where
fromInteger c = [fromInteger c]

negate fs = map negate fs

(f:ft) + (g:gt) = f+g : ft+gt
fs + [] = fs
[] + gs = gs

(f:ft) * gs@(g:gt) =
dropZeros \$ f*g : ft*gs + [f]*gt
_ * _ = []

instance Fractional Poly where
fromRational c = [fromRational c]

(0:ft) / gs@(0:gt) = ft/gt
(0:ft) / gs@(g:gt) = 0 : ft/gs
(f:ft) / gs@(g:gt) = f/g : (ft-[f/g]*gt)/gs
[] / (0:gt) = []/gt
[] / (g:gt) = []
_ / _ = error "improper polynomial division"

-- Polynomial composition
(f:ft) # gs@(0:gt) = f : gt*(ft#gs)
(f:ft) # gs@(g:gt) = [f] + gs*(ft#gs)
[] # _ = []
(f:_) # [] = [f]

-- Polynomial integration
integ fs = dropZeros \$
0 : zipWith (/) fs (countFrom 1)

-- Polynomial differentiation
diff (_:ft) = zipWith (*) ft (countFrom 1)
diff _ = [0]

countFrom n = n : countFrom (n+1)

dropZeros = foldr f [] where
f elem acc | elem == 0 && null acc = acc
| otherwise             = elem:acc

-- Code for root finding

newton :: Double -> (Double -> Double) -> [Double]
newton initial f = iterate go initial where
go x = x - (approx / df)
where approx = f x
df = (f (x + delta) - approx) / delta
delta = 0.001
pickValue :: Double -> (a -> Double) -> [a] -> a
pickValue tol f =
-- give up after 1000 tries
foldr go err . take 1000
where
go x x' = if abs (f x) <= tol then x else x'
err = error "can't converge"
findValue :: Double ->
(a -> Double) ->
(Double -> a) ->
a
findValue tol p fnc =
-- we use an initial guess of 1
pickValue tol p \$ map fnc \$ newton 1 (p . fnc)
```