module Data.Algorithm.CubicSpline(cubicSplineCoefficients) where
import Prelude hiding (tail, init, head, last, minimum, maximum, foldr1, foldl1, (!!), read)
import Control.Arrow
import Numeric.LinearAlgebra
import Data.List(unfoldr)
import Safe
type PolyCos = [Double]
cubicSplineCoefficients :: [(Double, Double)] -> [PolyCos]
cubicSplineCoefficients xs = chunkBy 4 . concat . toLists $ linearSolve matrix solution
where
(matrix,solution) = (fromLists *** trans . fromLists . (:[])) .
unzip . map (first fillOut) . (extraConditions ++) . initNote "cubicSpline" . initNote "cubicSpline" . concat .
zipWith (map . first . (++)) (iterate (++[0,0,0,0]) []) $
map genEquations (xs `zip` drop 1 xs)
genEquations :: ((Double, Double),(Double,Double)) -> [(PolyCos, Double)]
genEquations ((x,y),(x',y')) = [ (zipWith (^) (repeat x) ([0..3]::[Int]) , y)
,(zipWith (^) (repeat x') ([0..3]::[Int]) , y')
,(deriv1 ++ map negate deriv1 , 0)
,(deriv2 ++ map negate deriv2 , 0)]
where deriv1 = [0, 1, 2*x', 3*x'^(2::Int)]
deriv2 = [0, 0, 2 , 6*x']
extraConditions = [ ([0, 0, 2, 6 * fst (headNote "cubicSpline" xs)], 0)
, (replicate (fullLen 4) 0 ++ [0, 0, 2, 6 * fst (lastNote "cubicSpline" xs)], 0)]
fillOut = take fullLen . (++ repeat 0)
fullLen = 4 * (length xs 1)
chunkBy :: Int -> [t] -> [[t]]
chunkBy n = unfoldr go
where go [] = Nothing
go x = Just $ splitAt n x