module Numeric.LazySplines where
import Data.List
import Control.Arrow
class Sampleable a where
at :: a -> Double -> Double
type Poly = [Double]
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
liftS :: Double -> Spline
liftS x = [(maxDuration,[x])]
instance Sampleable PolySegment where
at (_,poly) pt = poly `at` pt
instance Sampleable Spline where
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) =
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'
| matchFirst = f totalDur dur $
match lastVal poly
| otherwise = matchScale lastVal dur $
f totalDur dur poly
match lastVal (_:xs) = lastVal:xs
match lastVal x = [lastVal]
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
doSplit = dur > minsize &&
abs (p t dur poly) > tol
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)
| dur <= 0 =
go (lastVal, time) (negate dur) fs
| otherwise =
(dur', poly') :
go (lastVal', time') chop' fs
where
poly = shiftBy chop oldPoly
dur = oldDur chop
chkPred d =
(abs $ p time d poly) < tol
dur' = lastDef dur .
takeWhile chkPred .
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 =
headDef poly .
map fst .
dropWhile chkPred $
zip polys (drop 1 polys)
where
polys = drop 2 . inits $ poly
chkPred (x,x') = p t dur x >
p t dur x'
headDef def [] = def
headDef _ (x:_) = x
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"
(f:ft) # gs@(0:gt) = f : gt*(ft#gs)
(f:ft) # gs@(g:gt) = [f] + gs*(ft#gs)
[] # _ = []
(f:_) # [] = [f]
integ fs = dropZeros $
0 : zipWith (/) fs (countFrom 1)
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
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 =
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 =
pickValue tol p $ map fnc $ newton 1 (p . fnc)