-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.CubicSpline
-- Copyright   :  (c) 2011 diagrams-lib team (see LICENSE)
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  diagrams-discuss@googlegroups.com
--
-- A /cubic spline/ is a smooth, connected sequence of cubic curves
-- passing through a given sequence of points.  This module implements
-- a straightforward spline generation algorithm based on solving
-- tridiagonal systems of linear equations.
--
-----------------------------------------------------------------------------
module Diagrams.CubicSpline.Internal
       (
         -- * Solving for spline coefficents
         solveCubicSplineDerivatives
       , solveCubicSplineDerivativesClosed
       , solveCubicSplineCoefficients
       ) where

import           Diagrams.Solve.Tridiagonal

import           Data.List

-- | Use the tri-diagonal solver with the appropriate parameters for an open cubic spline.
solveCubicSplineDerivatives :: Fractional a => [a] -> [a]
solveCubicSplineDerivatives :: [a] -> [a]
solveCubicSplineDerivatives (a
x:[a]
xs) = [a] -> [a] -> [a] -> [a] -> [a]
forall a. Fractional a => [a] -> [a] -> [a] -> [a] -> [a]
solveTriDiagonal [a]
as [a]
bs [a]
as [a]
ds
  where
    as :: [a]
as = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) a
1
    bs :: [a]
bs = a
2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) a
4 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
2]
    l :: Int
l  = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ds
    ds :: [a]
ds = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
f ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a] -> a
forall a. [a] -> a
last [a]
xs]) (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs)
    f :: a -> a -> a
f a
a a
b = a
3a -> a -> a
forall a. Num a => a -> a -> a
*(a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
b)

solveCubicSplineDerivatives [a]
_ = [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"argument to solveCubicSplineDerivatives must be nonempty"

-- | Use the cyclic-tri-diagonal solver with the appropriate parameters for a closed cubic spline.
solveCubicSplineDerivativesClosed :: Fractional a => [a] -> [a]
solveCubicSplineDerivativesClosed :: [a] -> [a]
solveCubicSplineDerivativesClosed [a]
xs = [a] -> [a] -> [a] -> [a] -> a -> a -> [a]
forall a. Fractional a => [a] -> [a] -> [a] -> [a] -> a -> a -> [a]
solveCyclicTriDiagonal [a]
as [a]
bs [a]
as [a]
ds a
1 a
1
  where
    as :: [a]
as = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) a
1
    bs :: [a]
bs = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
l a
4
    l :: Int
l  = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs
    xs' :: [a]
xs' = [a] -> [a]
forall a. [a] -> [a]
cycle [a]
xs
    ds :: [a]
ds = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
l ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
f (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
xs') (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [a]
xs')
    f :: a -> a -> a
f a
a a
b = a
3a -> a -> a
forall a. Num a => a -> a -> a
*(a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
b)

-- | Use the cyclic-tri-diagonal solver with the appropriate parameters for a closed cubic spline.
solveCubicSplineCoefficients :: Fractional a => Bool -> [a] -> [[a]]
solveCubicSplineCoefficients :: Bool -> [a] -> [[a]]
solveCubicSplineCoefficients Bool
closed [a]
xs =
    [ [a
x,a
d,a
3a -> a -> a
forall a. Num a => a -> a -> a
*(a
x1a -> a -> a
forall a. Num a => a -> a -> a
-a
x)a -> a -> a
forall a. Num a => a -> a -> a
-a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
da -> a -> a
forall a. Num a => a -> a -> a
-a
d1,a
2a -> a -> a
forall a. Num a => a -> a -> a
*(a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
x1)a -> a -> a
forall a. Num a => a -> a -> a
+a
da -> a -> a
forall a. Num a => a -> a -> a
+a
d1]
    | (a
x,a
x1,a
d,a
d1) <- [a] -> [a] -> [a] -> [a] -> [(a, a, a, a)]
forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 [a]
xs' ([a] -> [a]
forall a. [a] -> [a]
tail [a]
xs') [a]
ds' ([a] -> [a]
forall a. [a] -> [a]
tail [a]
ds')
    ]
  where
    ds :: [a]
ds | Bool
closed    = [a] -> [a]
forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivativesClosed [a]
xs
       | Bool
otherwise = [a] -> [a]
forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivatives [a]
xs
    close :: [a] -> [a]
close [a]
as | Bool
closed    = [a]
as [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a] -> a
forall a. [a] -> a
head [a]
as]
             | Bool
otherwise = [a]
as
    xs' :: [a]
xs' = [a] -> [a]
forall a. [a] -> [a]
close [a]
xs
    ds' :: [a]
ds' = [a] -> [a]
forall a. [a] -> [a]
close [a]
ds